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Chapter  1 
Introduction 


Array  processing  has  been  an  important  part  of  signal  processing  in  the  past  few 
years.  It  can  be  defined  as  the  use  of  an  array  of  sensors  to  receive  and  sample  signals.  It 
then  processes  the  received  data  in  order  to  obtain  some  information  about  the  sampled 
signals.  Some  of  that  information  includes  the  number  of  sources,  the  locations  of  those 
sources  and  the  type  of  signals  being  emitted.  The  array  consists  of  sensors  located  at 
different  points  in  space  with  respect  to  a  reference  point.  There  are  two  types  of  sensor 
array  systems:  active  and  passive.  In  the  active  case,  a  known  waveform  is  generated  and 
reflected  by  the  target.  In  the  passive  case,  the  signal  received  by  the  array  is  generated 
by  the  target. 

Direction  of  Arrival  (DOA)  denotes  the  direction  from  which  the  wave  fields 
arrive  at  the  sensor  array.  The  goal  in  DOA  detection  and  estimation  is  to  accurately 
determine  the  number  of  sources  producing  waveforms  and  the  locations  of  those 
sources.  Some  of  its  applications  include  cellular  communications,  air  traffic  control, 
seismology,  sonar  and  bioengineering. 

Many  algorithms  have  been  proposed  to  solve  the  DOA  estimation  problem.  The 
Beamforming  method  [1]  is  one  of  the  simplest  techniques  for  DOA  estimation.  This 
method  consists  in  computing  the  output  signal  using  a  weighted  sum  of  the  sensor 
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outputs.  The  DOA  is  determined  by  locating  the  peaks  in  the  power  spectrum  of  the  array 
output.  This  technique  has  poor  resolution  and  requires  a  great  number  of  sensors  for  a 
good  estimation  of  the  DOA. 

The  Maximum  Likelihood  approach  [2]  was  proposed  by  Bohme  in  1984.  Using 
this  method,  the  signals  are  modeled  as  unknown  deterministic  elements,  rather  than 
random  signals.  The  Maximum  Likelihood  Estimator  of  the  DOA  is  obtained  by 
searching  over  the  array  manifold  for  the  D  steering  vectors  that  form  a  signal  subspace 
that  is  closest  to  the  array  data  input  matrix.  Closeness  is  measured  by  the  modulus  of  the 
orthogonal  projection  of  the  input  vectors  onto  that  subspace.  This  technique  has  high 
resolution,  but  is  very  computationally  intensive. 

The  Multiple  Signal  Classification  algorithm  (MUSIC)  [3]  was  proposed  by 
Schmidt  in  1987.  It  is  one  of  the  most  popular  techniques  for  DOA  estimation  because  it 
has  a  high  resolution,  but  is  also  less  computationally  intensive  than  the  Maximum 
Likelihood  technique.  MUSIC  makes  use  of  the  eigenstructure  of  the  input  covariance 
matrix  by  separating  the  eigenvectors  into  orthogonal  subspaces.  The  MUSIC  approach 
was  first  proposed  for  the  narrowband  case  [4-6].  It  can  also  be  used  to  determine  the 
DOA  of  wideband  signals  after  separation  of  those  signals  into  narrowband  components 
[7-9], 

In  this  thesis,  a  parallel  architecture  for  DOA  estimation  of  wideband  sources  is 
devised.  The  algorithm  used  is  called  the  Coherent  Signal-Subspace  (CSS)  method  and 
was  proposed  by  Wang  and  Kaveh  [10].  This  technique  separates  the  wide  frequency 
band  into  narrowband  components.  The  covariance  matrices  for  all  the  frequency 
components  are  estimated  and  combined  to  form  a  single  focused  covariance  matrix.  The 
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narrowband  MUSIC  algorithm  can  then  be  applied  to  the  resulting  focused  covariance 
matrix. 

In  this  chapter,  the  narrowband  signal  model  is  given  and  the  MUSIC  algorithm 
is  presented.  The  wideband  data  model  is  also  presented  along  with  the  Coherent  Signal 
Subspace  method  for  wideband  sources. 

1.1  Narrowband  sources 

1.1.1  Signal  model 


Let  us  consider  an  array  of  M  identical  sensors.  Assuming  that  D  narrowband 
sources  (D  <  M)  with  identical  bandwidth  impinge  on  the  array  from  directions ...,0D , 

the  received  signal  at  the  ith  sensor  can  be  expressed  as: 


x.-(0"Sa*s*  (r-Tr)+n,(0 


m 


(i.i) 


where  sk  (t)  is  the  signal  at  a  reference  point,  xik  is  the  propagation  delay  between  the 

reference  point  and  the  ith  sensor  for  the  k"‘  source,  aik  is  the  impulse  response  of  the 

ith  sensor  to  the  kth  source,  and  ni  (f  )is  the  additive  noise  at  the  ith  sensor.  Since  sk  (f ) 

is  a  narrowband  process,  we  can  approximate  the  time  delay  by  a  phase  shift  [11], 
Equation  (1.1)  can  be  rewritten  as: 


(1.2) 


T 


ik 


.  sin0 
=  A - 


c 


(1.3) 
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where  A  is  the  distance  between  the  reference  point  and  the  sensor,  6  is  the  direction  of 
arrival  and  c  is  the  velocity  of  the  wave.  The  model  used  to  represent  the  output  of  the  M 
sensors  is: 


X(f)=A(0)S(f)+N(f)  (1.4) 

A(0)-[a(0,),  3(0,),  •••  a  (9d  )]  (1.5) 


a 


aM  (dk)e-JW°rM 


(1.6) 


Matrix  A  is  called  the  direction  matrix  and  each  of  the  column  vectors  are  the  direction 
vectors  of  the  sources.  Assuming  that  the  noise  is  independent  of  the  signal  with  zero 
mean,  the  covariance  matrix  can  be  expressed  as: 


R„-£[X(/)X"(r)] 

(1.7) 

R„  -  A (6)E [S«)S“ (r)]  A” (0)  +  E [N (»)N“  (»)] 

(1.8) 

If  the  output  array  vector  X(t)  is  observed  over  K  subintervals  of  duration  A T  seconds 
each,  the  covariance  matrix  [10]  can  be  expressed  as  the  snapshot  averaged  cross-product 
of  xt(f)  : 

R«  =^2xi(0xf  (0  (L9) 

a  m 

1.1.2  Multiple  signal  classification  algorithm: 

The  MUSIC  (Multiple  Signal  Classification)  algorithm  is  a  high  resolution 
technique  based  on  some  properties  of  the  input  covariance  matrix.  Following  the  signal 
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model  discussed,  the  received  data  vector  can  be  expressed  as  a  linear  combination  of  the 
incident  waveforms  and  noise  [12]: 

X(f)=A(0)S(f)+N(f)  (1.10) 


or 

X  =  AS  +  N  (1.11) 

where  s  =  [Sl  (r),s2  (r),...,sD  (r)]r  is  the  incident  vector  and  n  =  [^  (f),n2  (?),.. ,,nD  (f)]r  is  the 

noise  vector. 

The  input  covariance  matrix  can  be  expressed  as: 


R  xx=E[XXH] 

=  A E  [SSH  ]  Ah  +  E  [NN"  ] 

R„  =  ARSiAH+a„2I 

where  R55  is  the  signal  correlation  matrix  E  [SSH  j  and  o2n 

eigenvalues  of  the  input  covariance  matrix  R  u  are  the  values 
following  equation: 


(1.12) 

(1.13) 

(1.14) 

is  the  noise  variance.  The 
that  satisfy  the 


R„-V  -o 


(1.15) 


We  can  rewrite  the  previous  equation  as 


AR„  A"  +  o  \  -  AI 


n  i 


=  0 


ar„a"-(a,-ct„i)i|.o 

Therefore,  the  eigenvalues  of  AR  A"  are 


(1.16) 

(1.17) 


V;  =  A.  -  ol 


(1.18) 
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Matrix  A  is  composed  of  linearly  independent  direction  vectors  and  therefore  has  full 
column  rank  D .  The  signal  correlation  matrix  is  non  singular  as  long  as  the  incoming 
signals  are  not  highly  correlated.  As  a  result,  if  the  number  of  signals  D  is  less  than  the 
number  of  array  sensors  M,  the  M  x  M  matrix  ARSJAH  is  positive  semi  definite  with 

rankD .  This  implies  that  the  matrix  ARsvA//  had  D  non-zero  eigenvalues.  From  (1.18), 
we  see  that  M  -  D  eigenvalues  of  Ru  are  equal  to  the  noise  variance  o2n  .  After  sorting 
the  eigenvalues  of  Ru  so  that  A,  is  the  largest  value  and  XM  the  smallest,  we  have 


Ad+1,...,Am  =a„2  (1.19) 

The  eigenvector  associated  with  the  eigenvalue  A(  is  the  vector  q;  such  that 

(R„--Mk-°  (1.20) 

For  the  eigenvectors  associated  with  th eM -  D  eigenvalues,  we  have 

(R„-K,i)q,-o  (i.2i) 

(122> 

(AR„AH+CT,;i-CT„2l)q,-0  (1.23) 

AR„A“q,-0  (1.24) 


Since  A  has  full  rank  and  RJS  is  non  singular,  we  have 


A"q;  =  0 


'  a11  (di  )q,  ‘ 

0‘ 

)q, 

0 

(1.25) 


(1.26) 
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This  means  that  the  eigenvectors  associated  with  the  M  -  D  smallest  eigenvalues  are 
orthogonal  to  the  direction  vectors  making  up  A .  These  observations  form  the  basis  of 
the  MUSIC  algorithm.  The  analysis  shows  that  the  eigenvectors  of  the  covariance  matrix 
belong  to  two  orthogonal  subspaces  called  the  signal  and  noise  subspaces.  We  can 
estimate  the  direction  of  arrival  by  finding  direction  vectors  which  lie  in  the  signal 
subspace.  These  vectors  are  the  direction  vectors  that  are  orthogonal  to  the  noise 
subspace.  To  search  the  noise  subspace,  we  form  a  matrix  containing  the  noise 
eigenvectors: 

V.-[qM  -  q„]  (1.27) 


Since  the  direction  vectors  of  the  incoming  signals  are  orthogonal  to  the  noise  subspace 
eigenvectors,  we  can  say: 

(1.28) 

where  6  is  the  direction  of  arrival  of  a  signal  component. 

The  direction  of  arrival  can  then  be  estimated  by  finding  the  peaks  of  the  MUSIC 
spectrum  given  by: 


*>(«)- 


a"  (e)v„v”  m 


-,0e[O,2ji] 


(1.29) 


The  D  largest  peaks  in  the  spectrum  will  correspond  to  the  directions  of  arrival  of  the 
signals  impinging  on  the  sensor  array. 


The  MUSIC  algorithm  can  be  summarized  as  follow: 

1 .  Collect  input  samples  xk  ( t ) 

2.  Estimate  the  covariance  matrix  given  by: 
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(*K(0 

a  m 

3.  Compute  eigenvalues  and  eigenvectors  of  the  covariance  matrix. 

4.  Find  number  of  sources  D. 

5.  Find  the  DOA  estimates  by  finding  the  D  largest  peaks  of  the  MUSIC  spectrum 
given  by: 

P^=aH  JdjvXm 

1.1.3  Estimating  number  of  sources: 

The  MUSIC  algorithm  depends  on  the  parameter  D,  the  numbers  of  wave  fronts 
impinging  on  the  sensor  array,  for  estimating  the  DOA  of  multiple  targets.  The  Minimum 
Description  Length  (MDL)  was  used  in  order  to  achieve  this  objective  [13].  The  MDL  is 
specified  by: 


MDL(D)= 


-2  log 


M 

JI 


1  \ 
(m-d) 


M 


hi 


M 

k=D+\ 


+  -(2Af  -£>)log(A0 


\ 


/ 


(1.30) 


where  N  is  the  number  of  samples,  X  represents  the  eigenvalues  of  the  covariance  matrix 
and  M  is  the  number  of  sensors.  The  number  of  sources  is  determined  by  finding  a  value 
of  D  which  minimizes  the  MDL  criterion.  The  maximum  number  of  sources  that  can  be 
estimated  is  M  - 1 . 
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1.2  Wideband  sources 


1.2.1  Signal  model 

Let  us  consider  as  before  an  array  of  M  sensors.  Assuming  that  D  wideband 
sources  with  identical  bandwidth  B  impinge  on  the  array  from  directions eue2,...,eD,  we 

can  use  Equation  (1.1)  to  represent  the  signal  received  at  the  ith  sensor.  However,  the 
time  delay  cannot  be  approximated  as  a  phase  shift  in  the  wideband  case.  Assuming  that 
the  signals  are  observed  over  a  finite  interval  T,  we  can  represent  the  signal  x;  by  a 
Fourier  series  [11]: 

l+m 

no- 1  vk,  y*-'  d-31) 

n=l 

where  Xi  (wn )  are  the  Fourier  coefficients  given  by: 

T  /  2 

xt(wn)~^n  f  (1.32) 

1  -Til 

wn  =  —n,  n  =  l,...,l  +  m  (1.33) 

where  wl  is  the  lowest  frequency  and  wl+m  is  the  highest  frequency  included  in  the 

bandwidth  B.  Assuming  that  the  observation  time  T  is  much  greater  than  the  propagation 
delay  across  elements  of  the  array,  we  can  use  a  phase  shift  as  an  approximation  of  the 
time  delay  in  the  Fourier  domain. 

X‘  M=y,*ike-iWntlksk  (wn  )+n,  K  )  (1-34) 

The  model  used  to  represent  the  output  vector  is: 

X(w,)-A(w.)S(WJ+N(w„)  (1.35) 
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A(wn)=  [»OlK).  *62  M’  •••  aeo(>0] 


(1.36) 


a 


00= 


oi  ( ek  yiwt^ 
aM(6k)e-^ 


(1.37) 


As  a  result  of  the  Fourier  transform  applied  over  a  time  segment  AT ,  the  array 
output  vector  is  decomposed  into  non-overlapping  narrowband  components.  The 
covariance  matrix  for  component  wn  can  be  expressed  as: 


Rxxiwn)=  M 


R«(w«) =  A(^n  )E  [s(w„  )Sw  (wn  )1  XH  (wn)+E  [N  (wn  )Nh  (wn ) 


(1.38) 

(1.39) 


This  covariance  matrix  can  also  be  expressed  as: 

(Wn  )  ■ =  y.  *k  (Wn  K  M  (1 .40) 


i.2.2  Coherent  signal  subspace  method  for  wideband  sources 

It  is  possible  to  combine  the  signal  subspaces  of  different  frequencies  in  order  to 
generate  a  single  subspace  that  will  allow  us  to  determine  the  correct  number  of  sources 
and  directions  of  arrival  [10].  The  matrix  R  can  be  used  to  find  the  final  DOA  estimates. 
This  matrix  can  be  formed  by: 
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(1.41) 


r-2t(w,)r„(wj)t»(w;) 

7-1 

where  J  is  the  number  of  narrowband  components.  The  matrices  T  are  called 
transformation  matrices  and  can  be  expressed  as 


a2f,  (wo  )!<h p  ) 


aup  (W0  )laMfi  (Wl  ) 


(1.42) 


where  wa  is  the  central  frequency  of  bandwidth  B,  /3  is  the  initial  DOA  value, 
and  aif3  (mv  J  is  the  ith  element  of  the  direction  vector  (wj 


The  coherent  signal  subspace  method  for  computation  of  DOA  (wideband  sources)  as 
proposed  by  Wang  &  Kaveh  [10]  can  be  summarized  as  follow: 

1 .  Collect  data  samples  and  convert  the  samples  into  frequency  domain  using  FFT. 

2.  Estimate  the  covariance  matrix  for  each  frequency  component  given  by: 

(Wn  )  ■ =  Y  At  (Wn  K  (Wn  ) 

^  m 

3.  Compute  eigenvalues  and  eigenvectors  of  the  covariance  matrix. 

4.  Find  initial  estimates  of  direction  of  arrival  by  computing  the  MUSIC  spectrum 
given  by: 

5.  Compute  transformation  matrix  focusing  on  central  frequency. 
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6.  Compute  eigenvalues  and  eigenvectors  of  the  focus  matrix. 

7.  Find  number  of  sources  D. 

8.  Find  the  final  DOA  estimates  by  finding  the  D  largest  peaks  of  the  MUSIC 
spectrum. 

1.3  Thesis  Organization 

The  rest  of  the  thesis  is  organized  as  follows.  In  Chapter  2,  we  present  some  of  the 
methods  that  can  be  used  to  compute  the  eigenvalues  and  eigenvectors  of  the  covariance 
matrix.  Chapter  3  gives  the  simulations  details  for  a  single  DSP  implementation  of  the 
Coherent  Signal  Subspace  method  and  their  results.  In  Chapter  4,  a  parallel  Coherent 
Signal  Subspace  algorithm  is  presented  and  its  performance  is  measured.  Chapter  5 
presents  the  conclusion  and  future  research  topics. 
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Chapter  2 


Eigendecomposition  using  Householder  Matrices 


The  eigendecomposition  problem  is  a  very  important  part  of  the  DOA  estimation 
algorithm.  As  explained  in  the  previous  chapter,  finding  the  eigenvalues  and  eigenvectors 
of  the  covariance  matrix  is  needed  to  construct  the  signal  and  noise  subspaces  that  the 
MUSIC  algorithm  will  use.  The  Householder  and  QR  algorithms  [14-15]  can  be  used  to 
compute  the  eigenvalues  and  eigenvectors  of  the  symmetric  covariance  matrix.  The 
Householder  algorithm  is  used  to  reduce  the  bandwidth  of  the  covariance  matrix  by 
transforming  it  into  tridiagonal  form.  The  eigenvalues  and  eigenvectors  can  then  be 
computed  using  the  QR  algorithm.  The  computational  cost  needed  for  computing  the 
eigenvalues  and  eigenvectors  of  a  tridiagonal  matrix  will  be  much  smaller  than  that  of  the 
original  symmetric  matrix.  This  chapter  provides  a  brief  narrative  outlining  the 
mechanics  of  the  Householder  and  QR  algorithms. 

2.1  Householder  algorithm 

The  Householder  algorithm  [14-15]  reduces  an  n  x  n  symmetric  matrix  A  to 
tridiagonal  form  by  n-2  orthogonal  transformations.  Each  of  the  transformations  is  used 
to  eliminate  part  of  a  column  and  its  corresponding  row.  By  similarity,  the  tridiagonal 
matrix  B  and  the  original  symmetric  matrix  A  have  the  same  eigenvalues.  Furthermore,  it 
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can  be  proven  that  if  (x,A)  is  an  eigenpair  of  A  and  N  is  an  orthogonal  matrix,  then 
(l\Hx,  A)  is  an  eigenpair  ofN^RN.  As  a  result,  it  is  possible  to  find  the  eigenvectors  of 
matrix  A  using  orthogonal  transformations  on  the  eigenvectors  of  B  . 

An  n  x  n  matrix  H  is  called  a  Householder  matrix  if 

H  =  I-2wwr  (2.1) 


This  equation  can  be  rewritten  as 


H^- 


uii 


u 


(2.2) 


where  w  =  — and  llull  =  u2 . 
ul  11  1,2 


Note  that  H  is  symmetric  and  orthogonal. 

Let  us  consider  H  such  that 
Hx  =  ael 

where  a  is  a  constant,  ej  is  the  unit  vector  [l,0,...,0]  and  x  =  \xv...,xn  J 

It  is  possible  to  set  u  =  x  +  /3e {  and  determine  /3  such  that  Hx  doesn’t  have  any 

component  in  the  direction  ofx ,  as  specified  in  (2.3). 
x  +  /3el 


(2.3) 


w  =  • 


u 


U 


\\x  +  l3e 
Hx  =  x  -  2wwrx 

Hx  =  x  -  2 


Hl2 


'  x  +  /3el  ^ 

(x  +  fiex  )J 

x  +  Be .  „ 

1  1  2  ) 

x  +  ^  2 

■X 


(2.4) 

(2.5) 

(2.6) 
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Hx  =  x  -  2 


Hx  =  x  -  2 


(x  +  fte{ )  (' 

xrx  +  ftt\ 

ix) 

x  +  fte: 

2 

x  +  ftel 

2 

x(xTx  +  /?efx)+  /tej  (xrx  +  /3e[x) 


x  +  j Sej 


Since  fte\x  =  ftxx 


xix + i  ie  - 2x  (ihe + y  2Pe  1  (H2 + ) 

x^x^  +  2/3xj  +  ft2  —  2 |x|^  -2ftxx  ^-2fttx  ^|x|^  +  >3^  ) 

X  +  P*l  I 


(2.10) 


HjO^i  (lxl£+fo) 

x  +  ^ej  J 


(2.11) 


For  the  assumption  in  (2.3)  to  be  valid,  we  need 

P2  =\\xfi 

p-±H2 

u-Mml  x.y 


The  sign  of  ft  will  be  chosen  to  be  the  same  as  the  sign  ofv, . 
As  a  result,  the  vector  w  can  be  expressed  as: 


(2.12) 

(2.13) 

(2.14) 
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\  +^gn(x1)||x||2i 


[ 

From  (2. 1 1),  we  have 
Hx  =  -sign  (xj  )||x||2  e, 
and 

a  =  -^n(x1)||x||2 


(2.15) 


(2.16) 


(2.17) 


With  the  appropriate  choice  of  /3  ,  the  Householder  matrix  H  is  used  to  zero  out  part  of  a 
vector.  For  example,  if 

u  =  [0  0  •••  xk+sign(xk)  ||x||2  xk+l  •••  x„]  (2.18) 

then 

T 

Hx  =  |\  x2  ■■■  xk_x  -sign(xk  )||x||2  0  Oj  (2.19) 

H  has  the  form: 


A 

where  H  is  the  Householder  matrix  such  that 
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H 


(2.20) 


xt 


=  «e, 

where  a  is  a  constant. 

The  tridiagonalization  process  begins  by  computing  the  1st  Householder  matrix,  where  x 
is  chosen  to  be  the  lower  n- 1  elements  of  the  first  column  of  A.  As  a  result,  the  lower  n- 2 
elements  are  eliminated. 


HjA  = 


1  0  0  0  0 
0 


H 


an  &2\ 


a 
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a 


n\ 


*nl 


k 

0 

0 


In 


(2.21) 


(2.22) 


The  first  orthogonal  transformation  is: 
f  an  k  0  01 


A!  =  HAH  = 


k 

0 

0 

0 


(2.23) 


The  next  step  involves  the  computation  of  the  second  Householder  matrix  using  the  lower 
n- 2  elements  of  the  second  column. 
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h2 


'10  0  O' 
0  10  0 
0  0 

0  0  h2 
0  0 


(2.24) 


This  process  is  repeated  n- 2  times  in  order  to  obtain  the  tridiagonal  matrix. 

The  flowcharts  in  Figures  2.1  and  2.2  illustrate  and  summarize  the  tridiagonalization 
process. 

2.2  QR  factorization 

The  basic  idea  behind  the  QR  algorithm  [14-15]  states  that  any  real  matrix  A  can 
be  decomposed  in  the  form 

A  =  QR  (2.25) 

where  Q  is  an  orthogonal  matrix  and  R  is  an  upper  triangular  matrix. 

As  for  the  tridiagonalization  process,  it  is  possible  to  use  Householder  transformations  to 
eliminate  the  elements  under  the  main  diagonal. 

Then 

HMlHn_2  H2H1A  =  R  (2.26) 

Using  the  fact  that  the  Householder  matrices  are  symmetric,  it  follows  that 
A  =  H1H2---Hn.1R  (2.27) 

From  (2.27),  we  get 

Q  =  HjH2  •••Hn_1  (2.28) 
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Start 


Figure  2.1 :  Flowchart  of  function  House(x,k,n,w)  which  computes  the  householder  vector  w  used  to 
eliminate  the  lowest  k  components  of  a  vector  x  of  length  n 
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The  QR  method  can  be  used  to  find  the  eigenvalues  of  a  tridiagonal  matrix  A  using  the 


following  procedure: 

A*-Q*R4  (2.29) 

A4+1  =  (2.30) 

Since  Qk  is  orthogonal,  (2.30)  can  be  rewritten  as 

R.-lQjA,  (2.31) 

and 

A,.,-[QjAiQ1  (2.32) 


The  matrix  A  will  converge  to  a  diagonal  matrix  that  contains  the  eigenvalues  of  the 
tridiagonal  matrix  because  of  the  properties  of  similar  matrices.  After  m  iterations,  the 
product  of  the  orthogonal  transformations  Q  =  Q1Q2-  -Qm  contains  the  eigenvectors  of 

the  tridiagonal  matrix.  The  eigenvectors  of  the  original  matrix  can  be  computed  using  the 
following  orthogonal  transformations: 

X  =  H1H2  H„.2Q  (2.33) 

where  X  is  the  matrix  that  contains  the  eigenvalues  of  the  original  symmetric  matrix,  Q 
is  the  matrix  that  contains  the  eigenvalues  of  the  tridiagonal  matrix,  and 
hph2,  ,h„_2  are  the  Householder  matrices  computed  in  the  tridiagonal  process.  The 
flowcharts  in  Figures  2.3  and  2.4  show  the  implementation  of  the  QR  algorithm. 
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No 


Figure  2.3:  Flowchart  of  function  QR  (A,Q,R)  which  computes  the  QR  factorization  of  a  matrix  A 
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Figure  2.4:  Flowchart  of  function  QRShift  (A,n)  which  computes  the  eigenvalues  of  tridiagonal  matrix  A 
after  n  iterations 
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Chapter  3 


DSP  Implementation  of  CSS  algorithm 

The  CSS  algorithm  is  based  on  matrix  computations  and  orthogonal 
transformations,  which  are  computationally  intensive.  To  implement  this  algorithm  in 
real-time,  hardware  capable  of  executing  millions  of  operations  per  second  is  required. 
Some  of  the  architectures  available  for  that  purpose  are  DSPs,  FPGAs  and  ASICs.  A 
general  purpose  DSP  was  selected  as  an  appropriate  platform  for  implementation  because 
of  its  ease  of  programming.  Also,  the  DSP  is  the  best  suited  for  matrix  and  floating  points 
computations.  This  chapter  provides  some  information  about  the  DSP  used  and  gives  the 
simulations  details.  Experimental  results  are  also  presented. 


3.1  Digital  Signal  Processor 

3.1.1  Microcontroller 

The  DSP  used  in  our  implementation  is  a  DIOPSIS™  740  by  Atmel.  DIOPSIS™ 
740  (D740)  is  a  high  performance  dual-core  processing  platform  for  audio, 
communication  and  beam-forming  applications,  integrating  a  mAgic  DSP  and  an 
ARM7TDMI  RISC  MCU[16].  The  D740  is  optimally  suited  for  floating  point 
applications  complex  domain  computations. 
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The  ARM7TDMI  embedded  microcontroller  core  is  a  member  of  the  Advanced 


RISC  Machines  (ARM)  family  of  general-purpose  32-bit  microprocessors.  The  ARM 
architecture  is  based  on  Reduced  Instruction  Set  Computer  (RISC)  principles.  It  is 
equipped  with  peripherals  such  as  serial  interfaces  (SPI,  USART),  timers,  watchdog 
audio  A/D  and  D/A  converters.  The  ARM7  is  equipped  with  32  Kbytes  of  on-chip 
memory  and  runs  at  half  the  clock  speed  of  mAgic. 

The  mAgic  DSP  is  the  VLIW  numeric  processor  of  the  D740.  It  operates  on  IEEE 
754  40-bit  extended  precision  floating-point  and  32-bit  integer  numeric  format.  mAgic  is 
capable  of  delivering  1  Giga  Floating-Point  Operations  Per  Second  (GFLOPS)  at  a  clock 
rate  of  100  MHz.  The  main  components  of  the  DSP  subsystem  are  the  core  processor,  the 
on-chip  memories  and  the  interfaces  to  and  from  the  ARM  subsystem.  The  operators 
block,  the  register  fde,  the  address  generation  unit  and  the  program  decoding  and 
sequencing  unit  compose  the  core  processor.  Figure  3.1  shows  the  architecture  of  the 
mAgic  DSP. 

The  mAgic  DSP  has  four  on-chip  memory  blocks:  the  program  memory,  the  data 
memory,  the  data  buffer,  and  the  dual  ported  memory  shared  with  the  ARM  processor. 

An  external  memory  interface  multiplexes  the  data  accesses  and  the  program  accesses  to 
and  from  the  external  memory.  The  program  memory  stores  the  Very  Long  Instruction 
Word  (VLIW)  program  to  be  executed  by  mAgic.  It  is  8  KWords  x  128-bit  single  port 
memory.  The  mAgic  internal  data  memory  is  made  of  three  memory  pages  that  are 
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Figure  3.1  mAgic  DSP  Block  Diagram  (Courtesy  of  Atmel  Corporation)  [16] 


partitioned  into  left  and  right  memory  banks.  Each  bank  contains  2  KWords  x  40-bit,  thus 
a  total  of  12  KWords.  Each  data  memory  bank  is  a  dual  port  memory  that  allows  four 
simultaneous  accesses,  two  read  and  two  write.  The  buffer  memory  is  a  dual  port 
memory,  and  contains  2  KWords  x  40-bits  for  each  bank.  The  buffer  memory  is  a  dual 
port  memory.  One  port  is  connected  to  the  core  processor  and  the  second  port  is 
connected  to  the  external  memory  interface.  The  maximum  external  memory  size  of 
mAgic  is  16  MWords  for  each  bank  (equivalent  to  32  Mword  or  160  Mbytes;  24-bit 
address  bus).  A  Direct  Memory  Access  (DMA)  controller  manages  the  data  transfer 
between  the  external  memory  and  the  buffer  memory.  The  last  memory  block  in  the 
address  space  of  the  mAgic  DSP  is  the  shared  memory  (PARM)  between  mAgic  and  the 
ARM  processor.  It  is  a  dual  port  memory  containing  512  Words  by  40-  bit  for  both  the 
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left  and  the  right  bank  (total  IK  by  40-bit).  This  memory  can  be  used  to  efficiently 
transfer  data  between  the  two  processors. 


3.1.2  Software 

Multicore  Application  Development  Environment  (MADE)  is  an  Integrated 
Development  Environment  (IDE)  that  can  be  used  to  develop  D740  applications  [16].  It 
includes  the  C  compilers  for  both  ARM  and  mAgic  DSP  based  on  GNU  compiling  tools 
named  as  GCC.  It  has  a  high-level  mAgic  DSP  macro-assembler/optimizer,  an  editor  tool 
with  syntax  highlighting.  It  also  includes  a  unified  debugging  environment  that  can  be 
interfaced  with  a  cycle  accurate  simulator,  or  with  a  development  board  through  the  GNU 
debugger  tool  called  GDB.  The  magic  C  compiler  contains  a  DSP  library  composed  of 
over  220  functions  such  as  Fast  Fourier  Transform  and  HR  and  FIR  filter  creation. 

MADE  can  run  exclusively  in  two  modes: 

1.  Cycle  Accurate  Simulator  mode,  i.e.  attached  to  the  cycle  accurate  simulator. 

2.  Diopsis  target  mode,  i.e.  connected  to  a  board  through  a  serial  link. 


3.1.3  Evaluation  board 

The  JTST  board  [16]  is  low-cost,  stand-alone,  general-purpose  module  that 
provides  the  appropriate  resources  in  order  to  evaluate  D740  DSP  performances  in  a  wide 
range  of  applications.  The  JTST  board  provides  the  following  resources  to  D740: 

•  Memories:  mAgic  SSRAM,  ARM  FLASH  and  ARM  SRAM 
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•  Stereo  Audio  CODECs  (4  in  +  4  out) 

•  Serial  I/O: 

•  1  USB  2.0  Full  (12  Mbps) 

•  2  RS232/LVTTL  asynchronous/synchronous  serial  I/O  lines 

•  2  SPI  serial  I/O  lines 

•  Reset  Logic  (Power  ON,  Push  Button,  WDG) 

•  PLL  -  Clock  Logic 

•  Configuration  DIP  SWITCH  &  Status  7-segment  Display 

•  Voltage  Regulators  5V/3.3V  &  5V/1 ,8V 

•  Connectors  (USART,  SPI,  USB,  PIO,  AUDIO,  JTAG-ICE,  EXT  CLK,  PSU) 
Figure  3.2  shows  the  layout  of  the  JTST  board. 

3.2  Simulation 


In  order  to  demonstrate  the  performance  of  the  DOA  algorithm  for  wideband 
signals,  a  linear  array  of  sixteen  equally  spaced  omni  directional  sensors  was  used.  The 


spacing  between  sensors  is  - ,  where  c  is  the  velocity  of  propagation  and  f0  is  the 

2fo 


central  frequency.  We  will  attempt  to  estimate  DOA  0,  and  02  of  two  source  signals. 
The  signals  are  stationary  zero  mean  band  pass  white  Gaussian  processes  with  central 
frequency  fa  =  100Hz  and  bandwidth  B  =  40Hz .  The  array  noise  is  also  stationary  zero 

mean  band  pass  with  the  same  pass  band  as  the  signal  with  a  SNR  of  lOdB  at  each 
sensor. 
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Figure  3.2  JTST  Layout  [16] 


3.2.1  Sampling 

As  mentioned  above,  the  sources  signals  and  the  noise  are  random  processes  with 
a  bandwidth  of  40  Hz.  To  avoid  aliasing,  the  sampling  frequency  should  be  higher  than 


the  Nyquist  frequency,  which  is  twice  the  bandwidth.  The  sampling  frequency  is  chosen 


to  be  300  Hz.  The  signal  will  be  observed  over  a  period  of  Ta  seconds  and  T0  will  be 

divided  into  k  =  64  segments.  On  each  of  those  segments,  the  array  output  along  with 
corresponding  noise  will  be  decomposed  into  narrowband  components  using  a  fast 
Fourier  transform.  We  will  choose  a  number  of  64  sample  points  per  segment  in  order  to 
use  a  64  point  FFT.  The  total  number  of  samples  taken  by  each  sensor  will  be  4096. 
Simulation  data  is  similar  to  the  method  described  by  Wang  &  Kaveh  [10]. 

3.2.2  Data  generation 

Both  the  signal  and  noise  are  white  Gaussian  processes.  They  can  then  be 
generated  by  passing  a  set  of  random  numbers  through  a  band  pass  filter.  Random 
numbers  can  be  generated  using  a  linear  congruence  algorithm,  which  is  explained  in 
[17].  The  band  pass  filter  is  generated  by  using  the  function  IIR1  in  the  DSP  library  and 
specifying  the  band  pass  frequencies.  The  signal  can  be  represented  as 
s(t)  =  [.Sj  (tc ) . . .  Sj  (t63 ) . . .  sn  (t0 ) . . .  sn  (t63 ) . . .  s64  (t0 ) . . .  ^64  (t63 )]  (3-1) 

where  each  sample  is  a  random  number  passed  through  a  band  pass  filter. 

The  array  output  can  then  be  generated  by  using  the  wideband  signal  model  in 
(1.35).  The  time  domain  samples  will  be  transformed  into  frequency  domain  by  applying 
a  64  point  FFT  to  each  of  the  64  segments.  The  contents  of  both  the  time  sample  vectors 
and  the  frequency  sample  vectors  after  FFT  are  shown  below. 


X(0' 

'SnM' 

s„(h) 

FFT 

- > 

J»(Wl) 

Sn(t  63  ). 
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where  sn(tk_l)  is  the  kth  sample  of  the  nth  segment  and  sn(wn)  represents  the  component 
of  the  nth  segment  at  frequency  wn . 

The  value  of  the  sensor  impulse  response  a  is  1  for  omni  directional  vectors.  The 
signal  received  at  each  array  is  the  frequency  domain  data  scaled  by  a  factor  of  eJW"T , 
which  is  the  phase  shift  at  that  particular  sensor,  and  then  increased  by  a  value  n(wn), 

which  is  the  noise  value.  The  following  pseudo-code  shows  the  generation  for  two 
signals  arriving  at  an  array  composed  of  sixteen  sensors: 

For  i  =  1  to  16 

For  i  =  1  to  64 

For  i  =  1  to  64 

Generate  random  signal  values  Sx  and  s2 
Filter  values  of  sx  and  s2 

Generate  random  noise  values  nx  and  filter  values 

Compute  time  delay  for  each  frequency 

Apply  FFT  to  Sx  and  S2 

Multiply  sx  and  S2  by  time  delay 

Add  nx  to  sx  and  nx  to  s2 

Add  sx  and  s2  to  obtain  x 

Send  resulting  matrix  to  external  memory 

End 

End 

End 
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The  output  of  the  sensor  array  is  a  4096x16  matrix  expressed  as 


xl(0)  ...  *f(0)  -  *|6(0) 

xj  (63)  •••  x"  (63)  •••  x} 6  (63) 

x^4(0)  x6"4(0)  *£j(0) 

^64 (63)  ...  Ag4(63)  ...  xg4(63) 


(3.3) 


where  x'"  (k )  is  the  mth  output  sample  of  frequency  wk  for  the  nth  segment. 

Due  to  the  memory  limitation  of  the  DSP,  the  output  data  will  be  saved  in  the  external 
memory  as  sixty  four  64x16  matrices,  where  each  matrix  represents  the  output  for  each  of 
the  64  segments.  This  is  done  using  a  memory  transfer  from  local  to  external  memory. 


3.2.3  Covariance  matrix  computation 

The  covariance  matrix  for  each  frequency  component  can  be  computed  using 
Equation  (1.40).  The  main  issue  resides  in  forming  the  matrix  X  containing  all  samples 
of  one  frequency  components.  The  contents  of  the  external  memory  are  the  following: 
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*1  (0)  -  *|6(0) 
x\(63)  ...  x{6(63) 


'  44(0)  -  *64(0)' 

*64  (63)  ...  *64(63) 

As  we  can  see,  the  samples  related  to  each  frequency  are  spread  across  the  64  matrices 
and  need  to  be  put  in  the  same  matrix.  This  can  be  done  by  using  the  memory  transfer 
function  to  transfer  each  of  the  data  matrices  and  then  store  the  row  corresponding  to  the 
frequency  needed.  The  following  diagram  shows  that  computation  process  for  the 
frequency  component  wa . 


External  memory 


Local  memory 


Local  memory 


(A)  (B) 

x\(0)  -  x;6(0) 

x!(63)  ...  xl6(63) 
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Using  the  periodicity  and  symmetric  properties  of  the  FFT,  it  is  possible  to  have  all  the 


information  needed  by  selecting  frequencies  wa  to  w32 .  The  process  above  is  repeated 

for  each  of  the  33  covariance  matrices.  The  following  pseudo-code  shows  the  covariance 
matrices  computation  process: 

For  i  =  0  to  32 
For  j  =  0  to  63 

•  th 

Transfer  ]  data  matrix  to  local  memory  matrix  A 

Set  jth  row  of  matrix  B  equal  to  ith  of  matrix  A 

Compute  covariance  matrix  using  matrix  B 
Send  covariance  matrix  to  external  memory 

End 

End 


After  execution  of  this  code,  all  the  covariance  matrices  are  stored  in  external  memory. 
3.2.4  Eigendecomposition  of  covariance  matrix 

The  eigenvalues  and  eigenvectors  of  a  symmetric  matrix  can  be  computed  using 
the  Householder  and  QR  algorithms,  as  indicated  in  the  previous  chapter.  However,  these 
algorithms  cannot  be  directly  applied  to  the  covariance  matrix  as  it  is  a  Hermitian  matrix, 
i.e.  the  complex  analog  of  a  symmetric  matrix.  Solving  this  problem  requires  the 
conversion  of  the  Hermitian  matrix  into  a  real  symmetric  matrix.  It  should  be  noted  that 
Hermitian  matrices  have  real  eigenvalues. 

Let  C  =  A  +  /B  where  C  is  an  n  x  n  Hermitian  matrix,  A  and  B  are  n  x  n  real 
matrices.  The  eigenvalue  problem  for  a  complex  matrix  is 
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(A  +  zB)(u  +  zv)  =  A(u  +  zv) 


(3.4) 


where  A  is  the  real  eigenvalue  of  matrix  C,  and  u  +  zv  is  the  complex  eigenvector  of 
matrix  C.  This  n  x  n  complex  eigenvalue  problem  is  the  equivalent  to  the  2 n  x  2 n  real 
eigenvalue  problem. 


A 

-B 

. 

u 

=  A 

u' 

B 

A 

y 

V 

(3.5) 


The  2 n  x  2 n  matrix  is  only  symmetric  if  C  =  A  +  zB  is  Hermitian.  Solving  for  the 
eigenvalues  and  eigenvectors  of  the  2 n  x  2 n  matrix  using  the  Householder  and  QR 
methods  will  yield  two  2 n  x  2 n  matrices  containing  its  eigenvalues  and  eigenvectors. 
There  are  2 n  eigenvalues  in  the  form  A, ,  A, , . . . ,  A(1 ,  A„  where  A, , . . . ,  A(1  are  the  eigenvalues 
of  the  original  matrix  C.  The  eigenvectors  matrix  contains  the  eigenvectors  in  the  form 


u 


i  Ti 


h  «i 


n  n 


n  n 


(3.6) 


Corresponding  to  an  eigenvalue  A„ ,  the  vectors 


'-V  " 

n 

.V 

and 

n 

u« 

are  both  eigenvectors. 


Thus,  the  eigenvectors  pair  for  the  Hermitian  matrix  can  be  (un  +  i\ n )  or  (-v(j  +iun) 


The  following  pseudo-code  shows  the  computation  process  for  the  eigenvalues  of 
eigenvectors  of  the  covariance  matrices: 
for  i  =  0  to  32: 

•  th 

Transfer  l  covariance  matrix  to  local  memory  matrix  A 
Convert  A  into  real  symmetrical  matrix  B 
Find  eigenvalues  and  eigenvectors  of  matrix  B 
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Form  complex  eigenvectors  of  matrix  A 

Sort  eigenvalues  and  corresponding  eigenvectors  from  largest  to  smallest 
Transfer  eigenvalues  and  eigenvectors  to  external  memory 


End 


3.2.5  Power  spectrum 

The  power  spectrum  is  computed  using  Equation  (1.29).  As  explained  in  Section 
1.1.2,  the  column  vectors  forming  the  noise  matrix  are  the  eigenvectors  associated  with 
the  M-D  lowest  eigenvalues,  where  M  is  the  number  of  sensors  and  D  is  the  number  of 
sources.  The  power  spectrum  needs  to  be  computed  to  obtain  an  initial  estimate  for  the 
DOA.  The  number  of  sources  has  been  determined  using  the  MDL  algorithm.  The 
following  pseudo-code  shows  the  computation  process  for  the  power  spectrum: 

Input  value  d  for  the  number  of  sources 

Create  matrix  v  of  size  1 6  x  d 

for  i  =  0  to  90,  form  v  using  m-d  eigenvectors 
Form  direction  vector  A  using  angle  i 
Compute  power  spectrum 
Store  value  of  power  spectrum  in  array  Pm 

End 

Find  d  peaks  in  array  Pm 


The  value  of  D  used  is  1  for  the  initial  DOA  estimates. 
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3.2.6  Coherent  signal  subspace  processing 


As  explained  in  Section  1.2.2,  the  matrix  pencil  is  formed  using  the  covariance 
matrices  computed  in  Section  3.2.3.  The  covariance  matrices  must  be  transferred  back  to 
local  memory  for  the  computation  of  matrix  R .  The  process  is  included  in  the  following 
pseudo-code. 


Set  R  =  0 
Set  Rn=  0 
For  i  =  0  to  32 

Transfer  covariance  matrix  Cov( wf  )  to  local  memory 

Compute  transformation  matrix  T  using  central  frequency  and  initial  DOA  estimate 
Compute  T  •  Cov( wi  )-TH 
Add  result  to  R 

End 

The  eigenvalues  and  eigenvectors  of  the  matrix  pencil  will  be  computed  using  the 
Householder  and  QR  algorithms.  Thereafter,  they  are  used  to  determine  the  final  DOA 
estimates. 

3.3  Simulation  results 

The  Figures  3.3  to  3.6  show  the  results  of  simulations  for  estimation  of  DOA  of 
two  sources  located  in  the  interval  [0,  90].  For  each  case,  the  initial  and  final  DOA 
estimates  are  shown. 
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Case  1:  0j  =50  and02  =20 


Figure  3.3:  Plot  of  power  spectrum  for  initial  DOA  estimates  of  angles  qi  =50°  and  g2  =20° 


Figure  3.4:  Plot  of  power  spectrum  for  final  DOA  estimates  of  angles  gi  _  50°  and  <9,  _  20” 
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Case  2:  dx  =  30°  and02  =  40° 


Figure  3.5:  Plot  of  power  spectrum  for  initial  DOA  estimates  of  angles  qx  =  30°  and#2  =  40° 


Figure  3.6:  Plot  of  power  spectrum  for  final  DOA  estimates  of  angles  Qx  =30°  and  02  =  40° 
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The  results  indicate  that  DOA  detection  and  estimation  was  successful.  However,  a  high 
number  of  iterations  (40)  were  required  by  the  QR  algorithm  for  a  good  approximation  of 
the  eigenvalues  and  eigenvectors. 


40 


Chapter  4 


Parallel  Architecture  for  Coherent  Signal  Subspace 

Algorithm 


The  coherent  signal  subspace  method  for  wideband  DOA  algorithm  uses  the 
MUSIC  algorithm  which  is  a  high  resolution  method  capable  of  yielding  very  accurate 
DOA  estimates.  However,  this  algorithm  is  computationally  very  intensive.  The  overall 
computational  time  can  be  significantly  decreased  through  the  use  of  parallel  systems. 
The  following  approach  is  used  in  the  implementation  process: 

•  The  areas  that  can  be  improved  are  identified. 

•  The  algorithms  used  are  converted  into  computationally  efficient  modules. 

•  The  optimized  modules  are  divided  into  parallel  processes. 

In  this  chapter,  the  coherent  signal  subspace  algorithm  is  implemented  on  parallel 
processors  to  improve  its  performance.  Sections  of  the  algorithm  that  do  not  lend 
themselves  to  parallelization  will  be  implemented  on  a  single  DSP. 

4.1  Covariance  matrix  computation 

In  Section  3.2.3,  a  method  for  computing  the  covariance  matrices  using  a  single 
DSP  was  presented.  This  method  can  be  improved  upon  by  using  one  DSP  to  compute 
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each  of  the  33  covariance  matrices.  All  33  DSPs  share  the  same  external  memory.  Figure 
4.1  shows  its  parallel  architecture. 

The  computation  process  for  the  covariance  is  described  below. 

After  the  data  generation  step,  the  contents  of  the  external  memory  are  the 
following: 

xl(0)  ...  x\6(0) 

4(63)  ...  4 6 (63) 

44(0)  ...  44(0) 

4>4  (63)  ...  4(63) 

Processors  can  initiate  a  DMA  transfer  from  the  external  memory  to  their  local  memory. 
The  address  used  during  the  transfer  is  the  address  of  the  row  containing  the  samples  at  a 
particular  frequency  and  the  number  of  elements  transferred  is  the  number  of  elements 
contained  in  the  row.  The  DMA  transfer  for  the  first  matrix  is  shown  below: 


External  memory  Local  Memory 

4(0)  •••  x\6 (0)  1  -*  [4(0)  •••  46(0)]  DSP  0 

4(32)  •••  4 6 (32)  -4:(32)  46(32)]  DSP  32 
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Figure  4.1:  Parallel  architecture  for  coherent  signal  subspace  algorithm 
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After  a  DMA  transfer  is  complete,  another  transfer  will  be  initiated  for  a  different 
segment  matrix.  This  process  is  repeated  for  all  64  matrices  (segments)  in  external 
memory.  The  pseudo-code  executed  by  the  ith  processor  is  shown  as: 


For  m  =  0  to  63 

•  th 

Go  to  address  l  of  m  matrix  in  external  memory  and  transfer  1 6  elements  to  vector  v 
in  local  memory 

th 

Transfer  vector  v  to  M  row  of  vector  A 

End 

Compute  covariance  matrix  using  matrix  A 
Send  covariance  matrix  to  external  memory 


The  execution  time  will  be  the  same  for  each  processor.  This  execution  time  can 
be  further  reduced  by  computing  only  the  lower  triangular  part  of  the  covariance  matrix. 
Since  the  covariance  matrix  is  Hermitian,  the  information  contained  in  the  lower 
triangular  part  is  sufficient  to  form  the  entire  matrix.  Using  a  single  DSP  approach,  a 
total  of  64  x  33  DMA  transfers  of  256  elements  each  were  required.  The  parallelized 
process  only  requires  64  transfers  of  16  elements  each,  which  should  result  in  increased 
performance. 


4.2  Eigenvalues  and  Eigenvectors  computation 

The  eigenvalues  and  eigenvectors  computation  process  for  the  covariance  matrix 
is  based  on  two  algorithms,  the  Householder  and  the  QR  methods,  which  require  the  use 
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of  Householder  matrices.  Using  certain  properties  of  the  Householder  matrices,  it  is 
possible  to  parallelize  the  process.  Let  A  be  an  n  x  n  matrix  and  H  be  an  n  x  n 
Householder  matrix. 

If  B  =  HA  then 


B  =  (I-2wwr)A  (4.1) 

B  =  A-2wwrA  (4.2) 

If  we  let  b  ;  be  the  jth  column  vector  of  B  and  a .  be  the  jth  column  vector  of  A,  then 


b,=(CA), 


b  .  =  Ha . 


b  =  (I  -  2wwr  )a . 


b;  =  a  j  -  2ww  a . 


(4.3) 

(4.4) 

(4.5) 

(4.6) 


since  w  is  a  n  x  1  vector,  the  product  wray.  will  result  in  a  real  value  that  can  be  called 
a  .  Equation  (4.6)  can  then  be  expressed  as 

b  ;  =  a^  -  2aw  (4.7) 

It  is  thus  possible  to  compute  each  column  of  the  product  HA  separately.  For  a  matrix  A 
of  size  n  x  n,  n  processors  can  be  used  in  parallel  to  perform  this  computation.  This 
observation  will  allow  us  to  create  parallel  architectures  for  both  the  Householder  and  QR 
methods. 
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4.2.1  Householder  Method 


The  Householder  method  consists  of  transforming  a  symmetric  matrix  into 
tridiagonal  form  using  the  following  transformations 

B  =  Hn_2 . . .  H2Hj  AHjH2  •  •  •  Hm_2  (4.8) 

The  orthogonal  transformation  will  be  accumulated  in  a  matrix  Q  in  order  to  recover  the 
eigenvectors  of  the  original  matrix  A.  Equation  (4.8)  can  be  rewritten  as 
B  =  QAQ  (4.9) 

It  was  shown  that  the  product  HA  could  be  computed  using  n  parallel  processors  if  H  was 
a  Householder  matrix.  As  a  result,  H2HjA  and  Hn_2...H2H,A  can  also  be  computed 
with  n  processors.  If  we  let  C  =  QA ,  then  C  can  also  be  expressed  as 
C  =  H„_2...H2H1A  (4.10) 

Matrix  Q  is  expressed  as 

Q  =  Hm_2  ...H2HjI  (4.11) 

where  I  is  the  identity  matrix.  This  implies  that  C  and  Q  can  be  computed  using  the 
parallel  approach  explained  above.  The  tridiagonal  matrix  will  be  the  result  of  the  product 
CQ  .  Figure  4.1  shows  a  flowchart  of  the  parallel  Householder  method. 

After  execution  of  the  Householder  algorithm,  matrix  A  is  the  tridiagonal  matrix  and 
matrix  Q  contains  the  product  of  all  the  orthogonal  transformations. 
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Figure  4.1 :  Flowchart  of  parallel  Householder  method 


4.2.2  QR  method 


The  QR  method  is  based  of  the  use  of  the  following  orthogonal  transformations 
R  =  H„_j  ...H2HjA  (4.12) 

and 

Q  =  H„_1...H2H1  (4.13) 


Using  the  scheme  explained  in  the  previous  section,  we  can  compute  R  and  Q  using  n 
parallel  processors.  The  next  step  involves  the  computation  of  the  matrix  A  =  RQ  and  the 
computation  of  a  matrix  X  containing  the  product  of  the  orthogonal  transformations 
QiQ2"-Qm-where  m  is  the  number  of  iterations.  These  computations  can  be  done  using  a 

single  DSP.  This  process  will  start  over  until  the  number  of  iterations  required  for  a  good 
approximation  has  been  reached.  The  Figure  4.2  shows  a  flowchart  for  the  QR  method. 

After  execution  of  the  QR  algorithm,  the  matrix  A  contains  the  eigenvalues  of  the 
original  matrix  along  its  diagonal  and  matrix  X  contains  the  eigenvectors  of  the 
tridiagonal  matrix.  The  eigenvectors  of  the  original  matrix  can  be  obtained  be 
multiplying  the  matrix  Q  obtained  in  the  Householder  process  and  matrix  X. 
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Figure  4.2:  Flowchart  of  parallel  QR  method 
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4.2.3  Power  Spectrum 


The  power  spectrum  needs  to  be  computed  for  every  angle  between  0  and  90 
degrees.  Each  spectrum  value  can  be  computed  by  sending  the  matrix  containing  the 
eigenvectors,  the  number  of  sources  and  the  angle  of  arrival  to  a  specific  processor.  The 
DSP  can  then  compute  the  power  spectrum  and  send  the  results  back  to  external  memory. 
By  sending  3  angles  values  to  30  of  the  33  available  DSPs,  it  is  possible  to  reduce  the 
time  needed  to  compute  the  power  spectrum  by  a  factor  of  30. 


4.3  Simulations 

The  previous  simulation  created  for  single  DSP  was  used  again  for  the  case  of  the 
parallel  DSP  architecture.  This  was  done  to  measure  DOA  and  compare  their 
performances.  The  parallel  architecture  was  simulated  using  a  single  DSP  by  executing 
the  parallel  processes  sequentially.  However,  the  performance  would  be  measured  using 
the  longest  running  process  in  that  sequence  in  terms  of  number  of  clock  cycles.  The 
purpose  of  the  simulations  was  to  detect  and  estimate  the  DOA  of  two  sources  located  at 
20  and  60  using  40  iterations  for  the  QR  algorithm.  Figure  4.3  shows  the  simulations 
results  .It  can  be  seen  that  both  techniques  yielded  the  same  results. 

The  performance  analysis  showed  that  the  most  computationally  intensive  tasks  were  the 
QR  algorithm,  the  Householder  algorithm  and  the  power  spectrum  computation.  The 
following  table  shows  the  comparison  between  the  performance  of  the  single  DSP  and 
the  performance  of  the  parallel  architecture. 
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Final  DOA  Estimate 


Figure  4.3:  Plot  of  power  spectrum  for  final  DOA  estimate  of  angles  6X  =  20  and02  =  60  using  single 

and  parallel  DSP  approach 


Task 

Single  DSP 

Parallel  architecture 

Cycles 

Seconds 

Cycles 

Seconds 

Covariance  matrix 

1400000 

0.014 

42000 

0.00042 

Householder 

2600000 

0.026 

30000 

0.0003 

QR 

100000000 

1 

2100000 

0.021 

Power  spectrum 

1400000 

0.014 

47000 

0.00047 

Total 

210000000 

2.1 

5300000 

0.053 

Table  4.1  :  Performance  results  for  single  DSP  and  parallel  architecture 
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The  total  number  of  cycles  for  both  the  single  DSP  and  the  parallel  architecture  includes 
tasks  that  could  not  be  parallelized  and  are  not  listed  in  the  table.  The  total  execution  time 
for  the  parallel  architecture  is  0.05  seconds,  compared  to  2.1  seconds  for  the  single  DSP. 
This  is  due  to  the  fact  that  the  execution  time  for  QR  algorithm,  which  accounts  for 
approximately  95%  of  the  whole  process,  was  significantly  reduced  using  the  parallel 
architecture. 
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Chapter  5 


Conclusion 


5.1  Summary 

Direction  of  Arrival  detection  and  estimation  is  a  critical  part  of  array  processing. 
This  thesis  focused  on  developing  architecture  for  the  DOA  detection  and  estimation  for 
wideband  sources.  The  coherent  signal  subspace  algorithm  was  chosen  for  the 
implementation  because  of  its  high  resolution;  the  platform  selected  was  an  Atmel 
Diopsis740  DSP.  We  first  explained  the  difference  between  the  narrowband  and 
wideband  MUSIC  algorithms.  We  then  showed  that  the  coherent  signal  subspace 
algorithm  was  mainly  based  on  orthogonal  transformations  and  matrix  computations, 
which  had  high  computational  requirements.  A  parallel  architecture  capable  of  improving 
the  performance  of  the  coherent  signal  subspace  algorithm  was  proposed.  The 
performance  of  the  proposed  architecture  was  measured  and  compared  with  the  single 
DSP  implementation.  The  results  showed  that  the  parallel  architecture  yielded  the  same 
results,  while  providing  superior  performance. 
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5.2  Future  work 


One  of  the  limitations  of  the  proposed  architecture  was  the  use  of  static  matrices. 
This  implied  that  the  number  of  sensors  and  sources  was  known  in  advance.  It  also 
implied  that  the  system  would  not  be  able  to  detect  a  greater  number  of  sources  without 
major  modifications  in  the  source  code.  Using  dynamic  matrices  would  allow  the  system 
to  easily  adapt  to  the  number  of  sources  to  be  detected. 

A  Digital  Signal  Processor  was  chosen  for  the  implementation  because  of  its 
flexibility  and  ease  of  programming.  Other  platforms  such  as  ASICs  and  FPGAs  could  be 
used.  The  same  algorithm  could  be  implemented  on  those  platforms  so  that  the  best 
system  for  DOA  estimation  can  be  determined.  Some  of  the  criteria  for  comparison  could 
be  performance,  cost,  size  and  energy  consumption. 
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Chapter  7 


Appendix 

System  Software 

The  source  files  are  provided  on  a  compact  disk 
All  files  are  listed  below  along  with  a  brief  description: 

•  single. c  -  This  file  contains  the  single  DSP  source  code. 

•  parallel. c  -  This  file  contains  the  parallel  architecture  source  code. 

•  test.c-  This  file  contains  the  source  code  needed  by  the  ARM  processor 
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single.c  -  This  file  contains  the  single  DSP  source  code. 

#include  "math.h" 

#include  "magic. h" 
include  "DSPlib.h" 

#define  Pi  3.1415 
#define  nss  2 

#defme  velocity  300000000 
#define  omega  1885 
#defme  dist  500000 
#defme  ncc  16 
float  Pm[91]={23}; 
float  eigenvalues[16]={333}; 
float  eigvects[32][16]={14}; 
float  cova[32][32]={55}; 

_ complex _ float  Part[15]  ={5}; 

_ complex _ float  eigvects2[16][16]  ={5}; 

_ complex _ float  En[  16]  [15]  ={15}; 

_ complex _ float  nos[64][16]  BUFF={5}; 

_ complex float  time_delay[16][2]  BUFF  ={2,2, 4, 5, 3, 6, 9, 8, 7, 7, 4, 1,1, 1,1, 5}; 

float  buffi  1  [32] [32]  ={45}  ; 
float  buffi  [32] [32]  ={23}  ; 
float  w[32] ; 
float  rand(); 

//random  number  generator 

int  main() _ attribute _ ((halt, naked)); 

int  random_seed  =0  ; 
int  siggen(); 

//data  generation  function 

_ complex _ float  conj( _ complex _ float  a); 

void  covl(); 

//  Covariance  matrix  computation  parti 
void  cov2(); 

//Covariance  matrix  computation  part2 
void  fourier(); 

void  house(float  x[32],int  k);  //Householder  vector  function 

void  eigl(int  f); 
void  QR(); 

//QR  decomposition 
void  sym(); 

//Householder  method 

void  QRshift(int  f);  //  QR 

method 

int  count4; 

void  eigfunc(); 

void  bubbleSort(float  vals[],  int  array_size); 
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// 


void  music2(); 

Power  spectrum 
float  Q[32][32]={55}; 
float  R[32][32]  ={55}  ; 
float  H[32][32]={55}; 

_ complex _ float  W[32]  ={ 

#include  "coeff.dat" 


complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 


float  xntl[64][16]  XM; 
float  xnt2[64][16]  XM; 
float  xnt3[64][16]  XM; 
float  xnt4[64][16]  XM; 
float  xnt5[64][16]  XM; 
float  xnt6[64][16]  XM; 
float  xnt7[64][  16]  XM; 
float  xnt8[64][16]  XM; 
float  xnt9[64][16]  XM; 
float  xntl0[64][16]  XM; 
float  xnt  11  [64]  [16]  XM; 
float  xntl2[64][16]  XM; 
float  xnt  13  [64]  [16]  XM; 
float  xntl4[64][16]  XM; 
float  xntl5[64][16]  XM; 
float  xntl6[64][16]  XM; 
float  xntl7[64][16]  XM; 
float  xntl8[64][16]  XM; 
float  xntl9[64][16]  XM; 
float  xnt20[64][16]  XM; 
float  xnt21[64][16]  XM; 
float  xnt22[64][16]  XM; 
float  xnt23[64][16]  XM; 
float  xnt24[64][16]  XM; 
float  xnt25[64][16]  XM; 
float  xnt26[64][16]  XM; 
float  xnt27[64][  16]  XM; 
float  xnt28[64][16]  XM; 
float  xnt29[64][16]  XM; 
float  xnt30[64][16]  XM; 
float  xnt31[64][16]  XM; 
float  xnt32[64][16]  XM; 
float  xnt33[64][16]  XM; 
float  xnt34[64][  16]  XM; 
float  xnt35[64][16]  XM; 
float  xnt36[64][16]  XM; 


59 


complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 


float  xnt37[64][16]  XM; 
float  xnt38[64][16]  XM; 
float  xnt39[64][16]  XM; 
float  xnt40[64][16]  XM; 
float  xnt41[64][16]  XM; 
float  xnt42[64][16]  XM; 
float  xnt43[64][16]  XM; 
float  xnt44[64][16]  XM; 
float  xnt45[64][16]  XM; 
float  xnt46[64][16]  XM; 
float  xnt47[64][  16]  XM; 
float  xnt48[64][16]  XM; 
float  xnt49[64][16]  XM; 
float  xnt50[64][16]  XM; 
float  xnt51[64][16]  XM; 
float  xnt52[64][16]  XM; 
float  xnt53[64][16]  XM; 
float  xnt54[64][16]  XM; 
float  xnt55[64][16]  XM; 
float  xnt56[64][16]  XM; 
float  xnt57[64][  16]  XM; 
float  xnt58[64][16]  XM; 
float  xnt59[64][16]  XM; 
float  xnt60[64][16]  XM; 
float  xnt61[64][16]  XM; 
float  xnt62[64][16]  XM; 
float  xnt63[64][16]  XM; 
float  xnt64[64][16]  XM; 
float  bcov  1  [  1 6]  [  1 6]  XM; 
float  bcov2[16][16]  XM; 
float  bcov3[16][16]  XM; 
float  bcov4[16][16]  XM; 
float  bcov5[16][16]  XM; 
float  bcov6[16][16]  XM; 
float  bcov7[16][16]  XM; 
float  bcov8[16][16]  XM; 
float  bcov9[16][16]  XM; 
float  bcovl0[16][16]  XM; 
float  bcov  1 1  [  1 6]  [  1 6]  XM; 
float  bcovl2[16][16]  XM; 
float  bcov  1 3  [  1 6]  [  1 6]  XM; 
float  bcovl4[16][16]  XM; 
float  bcov  1 5  [  1 6]  [  1 6]  XM; 
float  bcovl6[16][16]  XM; 
float  bcovl7[16][16]  XM; 
float  bcov  1 8  [  1 6]  [  1 6]  XM; 
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_ complex _ float  bcovl9[16][16]  XM; 

_ complex _ float  bcov20[16][16]  XM; 

_ complex _ float  bcov21[16][16]  XM; 

_ complex _ float  bcov22[16][16]  XM; 

_ complex _ float  bcov23[16][16]  XM; 

_ complex _ float  bcov24[16][16]  XM; 

_ complex _ float  bcov25[16][16]  XM; 

_ complex _ float  bcov26[16][16]  XM; 

_ complex _ float  bcov27[16][16]  XM; 

_ complex _ float  bcov28[16][16]  XM; 

_ complex _ float  bcov29[16][16]  XM; 

_ complex _ float  bcov30[16][16]  XM; 

_ complex _ float  bcov3 1  [  1 6]  [  1 6]  XM; 

_ complex _ float  bcov32[16][16]  XM; 

_ complex _ float  bcov33[16][16]  XM; 

int  main() 

{ 

siggen(); 
fourier(); 
count4  =0; 

covl(); 

cov2(); 

eigl(0); 

sym(); 

QRshift(30); 

eigfunc(); 

bubbleSort(eigenvalues,  32); 

music2(); 

focus(); 

eig2(); 

sym(); 

QRshift(30); 

eigfunc(); 

bubbleSort(eigenvalues,  32); 

music2(); 

return  0; 

} 

_ complex _ float  conj( _ complex _ float  a){ 

_ complex _ float  b; 

_ real _ b= _ real _ a; 

_ imag _ b  =  -1* _ imag _ a; 

return  b;} 
float  rand() 

{ 

float  test; 
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random_seed  =  randomseed  *  25173  +13849; 

//  return  (unsigned  int)(random_seed  ); 

random_seed=(random_seed  %  1024); 
test  =  ((float)(random_seed))/1024; 
return  (test); 

} 

int  siggen()  { 

float  azi[2]={9,50}; 
float  sig_pow[2]={l,l}; 
float  nos_pow=0. 1 ; 
int  count, count2,count3; 
float  temp; 

int  index, index2,index3; 

_ complex _ float  sig_temp; 

_ complex _ float  sig[64][2]; 

nos_pow=sqrt(nos_pow) ; 

for  (count  =0;count<nss;count++){ 

azi[count]=azi[count]  *Pi/ 180; 

} 

nos_pow  =  sqrt(nos_pow); 

for  (count2  =0;count2<ncc;count2++){ 

for  (eount=0;count<nss;count++){ 

temp  =  sin(azi[count])*omega*dist*count2/velocity; 
//temp  =  temp*omega*count2*dist; 

_ real _ time_delay[count2][count]=  cos(temp); 

_ imag _ timedelay  [count2]  [count]=  sin(temp)*- 1 ; 

// 

} 

} 

//  index2=index*64; 

//  index3=  index2+64; 

for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<16;count2++){ 
nos[count]  [count2]=rand(); 
nos[count]  [count2]=nos[count]  [count2]  *nos_pow; 

} 

} 

for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<2;count2++){ 
sig[count]  [count2]=rand(); 

sig[count][count2]=sig[count][count2]*sig_pow[count2]; 

} 

} 

for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<16;count2++){ 
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sig_temp  =0; 

for  (count3  =  0;count3<2;count3++){ 
sig_temp 

=sig_temp+sig[count][count3]*time_delay[count2][count3]; 

} 

nos[count][count2]  =nos[count][count2]  +sig_temp; 

} 

} 


P32XM_B(nos,xntl ,  1 024); 
return  0; 

} 

void  fourier()  { 

int  countl,count2,count3; 

_ complex _ float  buffi  [64]; 

_ complex _ float  buff2[64]; 

_ complex _ float  tempi [64]; 

XM2P3_B(nos,xntl ,  1 024); 

for  (count  1  =0;countl<16;countl++){ 

for  (count2  =0;count2<64;count2++){ 
buffi  [count2]=nos[count2]  [count  1  ] ; 

//~//  testl[count2]  =buffl[count2]; 

} 

fft64(&W[0],  &buffl [0],  &templ[0],  &buff2[0]); 
for  (count2  =0;count2<64;count2++){ 
nos[count2]  [count  1  ]  =buff2  [count2] ; 

} 

} 

P32XM_B(nos,xntl ,  1 024); 

} 

void  covl(){ 

_ complex _ float  temp[64][16]; 

int  count, count2,count3; 

for  (count=0;count<64;count++){ 

XM2P3_B(nos,xntl ,  1 024); 

for  (count3  =0;count3<16;count3++){ 

temp[count][count3]=nos[count4][count3]; 

} 

} 

for  (count2  =0;count2<64;count2++){ 

for  (count  =0;count<16;count++){ 

nos[count2]  [count]=temp[count2]  [count] ; 

} 

} 
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} 

void  cov2()  { 

int  count, count2,count3; 

_ complex _ float  bcov[16][16]; 

for  (count  =0;count<16;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 
bcov[count]  [count2]=0; 

} 

} 

for  (count=0;count<64;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 

for(count3=0 ;  count3<  1 6 ;  count3++)  { 

bcov[count2][count3]=bcov[count2][count3]+(nos[count][count2]*conj(nos[coun 

t][count3])); 

} 

} 

} 

for  (count  =0;count<16;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 

bcov[count][count2]=bcov[count][count2]/64; 
nos[count]  [count2]=bcov[count]  [count2] ; 

} 

} 

P32XM_B(nos,bcovl,64); 

} 

void  eigl(int  f){ 
int  x,y; 
switch(f)  { 

case  0:  XM2P3_B(nos,bcovl,64); 
break; 

case  1 :  XM2P3_B(nos,bcov2,64); 
break; 

case  2:  XM2P3_B(nos,bcov3,64); 
break; 

case  3:  XM2P3_B(nos,bcov4,64); 
break; 

case  4:  XM2P3_B(nos,bcov5,64); 
break; 

case  5:  XM2P3_B(nos,bcov6,64); 
break; 

case  6:  XM2P3_B(nos,bcov7,64); 
break; 

case  7:  XM2P3_B(nos,bcov8,64); 
break; 

case  8:  XM2P3_B(nos,bcov9,64); 
break; 
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case  9:  XM2P3_B(nos,bcov  10,64); 


break; 

case  10: 

XM2P3_B(nos,bcovl  1,64); 
break; 

case  1 1 : 

XM2P3_B(nos,bcov  12,64); 
break; 

case  12: 

XM2P3_B(nos,bcov  1 3 ,64); 
break; 

case  13: 

XM2P3_B(nos,bcov  1 4,64); 
break; 

case  14: 

XM2P3_B(nos,bcov  1 5,64); 
break; 

case  15: 

XM2P3_B(nos,bcov  16,64); 
break; 

case  16: 

XM2P3_B(nos,bcov  17,64); 
break; 

case  17: 

XM2P3_B(nos,bcovl8,64); 

break; 

case  18: 

XM2P3_B(nos,bcov  19,64); 
break; 

case  19: 

XM2P3_B(nos,bcov20,64); 

break; 

case  20: 

XM2P3_B(nos,bcov2 1 ,64); 
break; 

case  21: 

XM2P3_B(nos,bcov22,64); 

break; 

case  22: 

XM2P3_B(nos,bcov23,64); 

break; 

case  23: 

XM2P3_B(nos,bcov24,64); 

break; 

case  24: 

XM2P3_B(nos,bcov25,64); 

break; 

case  25: 

XM2P3_B(nos,bcov26,64); 

break; 

case  26: 

XM2P3_B(nos,bcov27,64); 

break; 

case  27: 

XM2P3_B(nos,bcov28,64); 

break; 

case  28: 

XM2P3_B(nos,bcov29,64); 

break; 

case  29: 

XM2P3_B(nos,bcov30,64); 

break; 

case  30: 

XM2P3_B(nos,bcov3 1 ,64); 
break; 

case  3 1 : 

XM2P3_B(nos,bcov32,64); 

break; 
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case  32:  XM2P3_B(nos,bcov33,64); 

break; 

} 

//-  //XM2P3_B(m,bcov  1,16); 
for  (x  =0;x<16;x++){ 

for  (y=0;y<16;y++){ 

cova[x][y]  = _ real _ nos[x][y]; 

} 

} 

for  (x  =16;x<32;x++){ 

for  (y=0;y<16;y++){ 

cova[x][y]  =( imag nos[x-8][y])*(-l); 

} 

} 

for  (x  =0;x<16;x++){ 

for  (y=16;y<32;y++){ 

cova[x][y]  =( imag nos[x][y-8]); 

} 

} 

for  (x  =16;x<32;x++){ 

for  (y=16;y<32;y++){ 

cova[x][y]  = _ real _ nos[x-8][y-8]; 

} 

} 

} 

void  house(float  x[],int  k){ 
float  s, gamma; 
int  i; 

gamma  =  0; 

for  (i  =  0;i<(k-l);i++){ 

w[i]  =  0; 

} 

for  (i  =(k-l);i<32;i++){ 
gamma  +=  x[i]*x[i]; 

} 

gamma  =  sqrt(gamma); 
s  =  x[k-l]*x[k-l]; 
s  =  sqrt(s); 
s  =  gamma+s; 
s  =  2*gamma*s; 
s  =  sqrt(s); 
if  (x[k-l]>=  0){ 
w[k-l]  =  x[k-l]+gamma; 
w[k-l]  =  w[k-l]/s; 

} 

else  { 
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w[k-l]  =  x[k-l]-gamma; 
w[k-l]  =  w[k-l]/s; 

} 

for  (i  =  k;i<32;i++){ 
w[i]  =  x[i]/s; 

} 

} 

void  QR()  { 
float  buff2[32]; 

_ complex _ float  buff[32][32]; 

//float  buff3[32][32]; 
int  i,j,x,y,z; 
for  (i  =  0;i<32;i++){ 

for  (j  =  0;j<32;j++){ 

if  (i  ==j) imag buflf[i][j]  =1; 

else _ imag _ buff[i][j]  =  0; 

Q[i][j]  = imag buff[i]  [j] ; 

} 

} 

for  (i  =  0;i<31;i++){ 

for  (j  =0;j<32;j++){ 
buff2[j]  =  cova[j][i]; 

} 

house(buff2,(i+ 1 )) ; 

//gamma  =  0; 

for  (  x  =0;x<32;x++)  { 

for(y  =  0;y<32;y++){ 

H[x][y]  =  w[x]*w[y]; 

H[x][y]  =  2*H[x][y]; 

H[x][y]  = _ imag _ buff[x][y]  -  H[x][y]; 

_ real _ buff[x][y]  =0; 

} 

} 

for  (  x  =0;x<32;x++)  { 

for  (  y  =0;  y<32;  y++)  { 

for  (  z  =0;  z<32;  z++)  { 

_ real _ buff[x][y]  +=  Q[x][z]*H[z][y]; 

}  "  . ' 

} 

} 

for  (  x  =0;x<32;x++)  { 

for(y  =  0;y<32;y++){ 

Q[x][y]  = _ real _ buff[x][y]; 

_ real _ buff[x][y]  =  0; 

} 
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for  (  x  =0;x<32;x++)  { 

for  (  y  =0;  y<32;  y++)  { 

for  (  z  =0;  z<32;  z++)  { 

_ real _ buff[x][y]  +=  H[x][z]*cova[z][y]; 


} 

for(  x  =0;x<32;x++){ 

for(y  =  0;y<32;y++){ 

cova[x][y]  = _ real _ buff[x][y]; 

R[x][y]  = _ real _ buff[x][y]; 


void  sym()  { 

_ complex _ float  buff[32][32]; 

float  buff2[32]; 
float  buff3[32][32]; 

//~  float  buff4[n2][n2]; 

//~  float  buff6[n2][n2]; 

//~  float  buff7[n2][n2]; 
float  y[32]; 
float  z[32]; 
float  delta; 
int  x,s,i,j,v; 
for  (  x  =0;x<32;x++)  { 

for  (  s  =0;s<32;s++){ 

if  (x  ==s) _ real _ buff[x][s]  =1; 

else _ real _ buff[x][s]  =0; 

buff5[x][s]  = _ real _ buff[x][s]; 

} 

} 

for  ( i  =0;i<14;i++){ 

f°r  ( j  =0;j<32;j++){ 

buff2[j]  =  cova[j][i]; 
zD']  =  0; 

y[j]  =  0; 

} 

delta  =0; 

house(buff2,(i+2)); 
for  (  x  =0;x<32;x++)  { 

for  (  s  =0;s<32;s++){ 

if  (x  ==s) _ real _ buff[x][s]  =1; 

else _ real _ buff[x][s]  =0; 
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_ imag _ buff[x][s]  =  w[x]*w[s]; 

_ imag _ buff[x][s]  =  2* _ imag _ buff[x][s]; 

_ imag _ buff[x][s]  = _ real _ buff[x][s]  - _ imag _ buff[x][s]; 

}  ~  ' . “ 

} 

for  (x  =  0;x<32;x++)  { 

for  (  s  =0;s<32;s++){ 
for  (v  =0;  v<32;  v++){ 

buff3  [x]  [s]+= _ imag _ buff[x]  [v]  *cova[v]  [s] ; 

} 

} 

} 

for  (x  =  0;x<32;x++)  { 

for  (  s  =0;s<32;s++){ 
for  (v  =0;  v<32;  v++){ 

cova[x][s]+=  buff3[x][v]* _ imag _ buff[v][s]; 

}  .  . 

: 

} 

} 

} 

void  QRshift(int  f)  { 

_ complex _ float  buff[32]  [32] ; 

float  buff7[32][32]; 
float  buff2[32]; 
int  x,y,i,z; 

for  (x  =0;x<32;x++){ 

for  (y  =0;y<32;y++){ 

_ imag _ buff[x][y]  =0; 

if  (x  ==y) _ real _ buff[x][y]  =1; 

else _ real _ buff[x][y]  =0; 

} 

} 

for  (i=0;i<f;i++){ 

QR(); 

for  (x  =0;x<32;x++){ 

for  (y  =0;  y<32;  y++)  { 
buff7[x][y]  =0; 

_ imag _ buff[x][y]  =  0; 

for  (z  =0;  z<32;  z++)  { 

buff7[x][y]  +=  R[x][z]*Q[z][y]; 

_ imag buff[x][y]+= real buff[x]  [z]  *  Q  [z]  [y  ] ; 

} 

//  if  (fabs(buff[x]  [y])<  1  e-6)  buff[x]  [y]=0; 

} 

} 


69 


for(  x  =0;x<32;x++){ 

for  (y  =0;  y<32;  y++)  { 

cova[x][y]  =buff7[x][y]; 

_ real _ buff[x][y]= _ imag _ buff[x][y]; 

//  buff6[x][y]= _ imag _ buff[x][y]; 

buff2 1  [x]  [y]= _ imag _ buff[x]  [y] ; 

} 

} 

} 

intj; 

for  (  x  =0;x<32;x+=2){ 
j=x/2; 

eigenvalues  [j  ]  =buff7  [x]  [x] ; 

} 

} 

void  eigfunc()  { 
int  x,s,j,y,z; 

float  test2; 
float  buff23[32][32]; 
for  (  x  =0;x<32;x++)  { 

for  (  s  =0;  s<32;  s++)  { 
buff23[x][s]  =  0; 
for  (  z  =0;  z<32;  z++)  { 

buff23  [x]  [s]+=  buff5  [x]  [z]  *buff2 1  [z]  [s] ; 

} 

//  float  test2  =fabs(buff23[x][s]); 

//  if  (fabs(test2)<le-4)  { 

//  buff23[x][s]  =0;} 

} 

} 

for(  x  =0;x<32;x++){ 

for(  y  =  0;y<32;y+=2){ 
j=y/2; 

eigv  ects  [x]  [j  ]  =buff23  [x]  [y] ; 

} 

} 

} 

void  bubbleSort(float  vals[],  int  array_size) 

{ 

int  i,  j,  k; 
float  tempi; 
float  temp2[32]; 

for  (i  =  (array_size  -  1);  i  >=  0;  i— ) 

{ 

for  (j  =  1;  j  <=  i;  j++) 
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{ 

if  (valsQ-1]  <  valsQ]) 

{ 

tempi  =  vals[j-l]; 
vals[j-l]  =  vals[j]; 
vals[j]  =  tempi; 

for(k  =0;k<32;k++)  { 
temp2  [k] =eigvects  [k]  [j  - 1  ] ; 

} 

for(k  =0;k<32;k++)  { 
eigvects[k]  [j- 1  ]=eigvects[k]  [j] ; 

} 

for(k  =0;k<32;k++)  { 
eigvects  [k]  [j  ]  =temp2  [k] ; 

} 

} 

} 

} 

for  (i=0;i<16;i++){ 

for  (j=0;j<16;j++){ 

_ real eigvects2  [i]  [j  ]  =  eigvects  [i]  [j  ] ; 

_ imag _ eigvects2[i][j]=  eigvects[i+8][j]; 

} 

} 

} 

void  music2()  { 

int  k,  count; 
int  l,v; 

float  temp,temp2; 
int  kk; 

_ complex _ float  sr; 

_ complex _ float  array[16]={6}; 

_ complex _ float  array2[16]={7}; 

// _ complex _ float  Part[15]  ={5}; 

_ complex _ float  Part2[15]  ={5} ; 

float  xxx[91]=  {5}; 

_ complex _ float  En2[15][16]  ={5} ; 

_ vector _ int  error; 

for  (1  =  0;l<16;l++){ 

for  (v  =0;v<7;v++){ 

En[l]  [v] =eigvects2  [1]  [v+ 1  ] ; 

En2  [v]  [1] =eigvects2  [1]  [  v+ 1  ] ; 

} 

} 

for  (k=0;k<=90;k++){ 
xxx[k]=k*Pi; 
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xxx[k]  =xxx[k]/180; 
temp  =  sin(xxx[k]); 
temp  =  temp*Pi; 

for  (count  =  0;count<16;count++){ 
temp2  =  temp*count; 

_ real _ array[count]=  cos(temp2); 

_ real _ array2[count]=  cos(temp2); 

_ imag _ array  [count]=  sin(temp2)*-l; 

_ imag _ array2[count]=  sin(temp2)*-l; 

} 

int  s,j; 

for  (s=0;s<7;s++){ 

Part2[s]  =0; 

for  (kk=0;kk<16;kk++){ 

Part2  [s] +=conj  (En2  [s]  [kk] )  *  array  [kk] ; 

} 

} 

for  <j =0  ;j  <7  ;j ++)  { 

Part[j]  =0; 

for  (kk=0;kk<16;kk++){ 

Part[j  ]+=conj  (array  2  [kk])  *En[kk]  [j] ; 

} 

} 

Pm[k]  =0; 
sr  =0; 

for  (kk=0;kk<7;kk++)  { 

sr+=Part[kk]  *Part2[kk] ; 

} 

Pm[k]  =  ( _ real _ sr  * _ real _ sr)+( _ imag _ sr* _ imag _ sr); 

Pm[k]=sqrt(Pm[k]); 

} 

} 

void  focus()  { 
float  fO; 

float  temp,temp2,temp3; 
int  t,i,j,i2,j2,k2,12,m2; 

//bl=10.5; 

for  (i  =0;i<nfib;i++)  { 
fbin[i]  =(i+l)*300/16; 

} 

for  (i  =0;i<ncc;i++){ 

for  (j=0;j<ncc;j++){ 

m[i]D]  =0; 
signal[i][j]  =0; 

} 

} 
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t  =  (nfib+l)/2; 
fO  =fbin[t-l]; 


for(i=0;i<nfib;i++)  { 

for(i2=0;i2<nfib;i2++)  { 

for(j  2=0  ;j  2<nfib  ;j  2++)  { 
tl[i2][j2]  =0; 
tPl[i2][j2]  =0; 
tcl[i2][j2]  =0; 
tc2[i2][j2]  =0; 
ttl[i2][j2]  =0; 

} 

} 

omg  =2*Pi*(fO-fbin[i]); 

bl=20; 

bl  =bl*Pi; 

bl=bl/180; 

temp  =sin(bl); 

temp  =sin(bl)/300000000; 

temp  =temp*500000; 

temp  =temp*omg; 

// 

for(k2=0;k2<ncc;k2++)  { 
temp2=  temp*k2; 

_ real _ tl[k2][k2]=  cos(temp2); 

_ imag _ tl  [k2][k2]=  sin(temp2)*- 1 ; 

} 

switch  (i)  { 

case  0: 

XM2P3_B(tpp,bcovl  ,256); 
break; 

case  1: 

XM2P3_B(tpp,bcov2,256); 

break; 

case  2: 

XM2P3_B(tpp,bcov3 ,256); 
break; 

case  3: 

XM2P3_B(tpp,bcov4,256); 

break; 

case  4: 

XM2P3_B(tpp,bcov5 ,256); 
break; 

case  5: 

XM2P3_B(tpp,bcov6,256); 
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break; 

case  6: 

XM2P3_B(tpp,bcov7,256); 

break; 

case  7: 

XM2P3_B(tpp,bcov8 ,256); 
break; 

case  8: 

XM2P3_B(tpp,bcov9,256); 

break; 

case  9: 

XM2P3_B(tpp,bcov  10,256); 
break; 
case  10: 

XM2P3_B(tpp,bcovl  1,256); 
break; 
case  1 1 : 

XM2P3_B(tpp,bcov  12,256); 
break; 
case  12: 

XM2P3_B(tpp,bcov  13,256); 
break; 
case  13: 

XM2P3_B(tpp,bcov  14,256); 
break; 
case  14: 

XM2P3_B(tpp,bcov  15,256); 
break; 
case  15: 

XM2P3_B(tpp,bcov  16,256); 
break; 
case  16: 

XM2P3_B(tpp,bcov  17,256); 
break; 
case  17: 

XM2P3_B(tpp,bcovl  8,256); 
break; 
case  18: 

XM2P3_B(tpp,bcov  19,256); 
break; 
case  19: 

XM2P3_B(tpp,bcov20,256); 
break; 
case  20: 

XM2P3_B(tpp,bcov2 1 ,256); 
break; 
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case  21: 

XM2P3_B(tpp,bcov22,256); 
break; 
case  22: 

XM2P3_B(tpp,bcov23 ,256); 

break; 

case  23: 

XM2P3_B(tpp,bcov24,256); 

break; 

case  24: 

XM2P3_B(tpp,bcov25 ,256); 

break; 

case  25: 

XM2P3_B(tpp,bcov26,256); 

break; 

case  26: 

XM2P3_B(tpp,bcov27,256); 

break; 

case  27: 

XM2P3_B(tpp,bcov28,256); 

break; 

case  28: 

XM2P3_B(tpp,bcov29,256); 

break; 

case  29: 

XM2P3_B(tpp,bcov30,256); 

break; 

case  30: 

XM2P3_B(tpp,bcov3 1 ,256); 

break; 

case  3 1 : 

XM2P3_B(tpp,bcov32,256); 

break; 

case  32: 

XM2P3_B(tpp,bcov3  3,256); 
break; 

} 

for(i2=0;i2<ncc;i2++)  { 

tPl[i2][i2]  =conj(tl[i2][i2]); 
tt  1  [12]  [i2]=t  1  [i2]  [12]  *tp  1  [12]  [12] ; 

} 

for(  i2  =0;i2<ncc;i2++){ 
for  (  j2  =0;  j2<ncc;  j2++)  { 

tcl[i2][)2]  =  0; 
for  (k2  =0;  k2<ncc;  k2++){ 

tc  1  [12]  [j2]+=  1 1  [i2]  [k2]  *tpp[k2]  [j2] ; 
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} 

} 

} 

for(  i2  =0;i2<ncc;i2++){ 
for  (  j2  =0;  j2<ncc;  j2++)  { 

tc2[i2][j2]  =  0; 
for  (k 2  =0;  k2<ncc;  k2++){ 

tc2[i2]  [j2]+=  tc  1  [i2]  [k2]  *tp  1  [k2]  [j2] ; 

} 

signal[i2][j2]  =  signal[i2][j2]+tc2[i2][j2]; 

} 

} 

} 

for(  i2  =0;i2<ncc;i2++){ 
signal[i2][i2]= _ real _ (signal[i2][i2]); 

} 

} 

void  eig2()  { 

int  i2,j2; 

for  (i2  =0;i2<16;i2++){ 
for  (J2=0;j2<16;j2++){ 

A[i2][j2]  =_real_  signal  [i2][j2]*1e8; 

} 

} 

for  (i2  =16;i2<32;i2++){ 

for  (J2=0;j2<16;j2++){ 

A[i2][j2]  = _ imag _ signal  [i2-l  6]  [j2]*le8; 

} 

} 

for  (i2  =0 ;  i2<  1 6 ;  i2++)  { 

for  (J2=16;j2<32;j2++){ 

A[i2][j2]  =( _ imag _ signal[i2]  ^2- 16])*(- 1)*  1  e8; 

} 

} 

for  (i2  =16;i2<32;i2++){ 

for  (J2=16;j2<32;j2++){ 

A[i2][j2]  = _ real _ signal[i2-16]|j'2-16]*le8; 

} 

} 

} 
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parallel.c  -  This  file  contains  the  parallel  architecture  source  code. 

#include  "math.h" 

#include  "magic. h" 

#include  "DSPlib.h" 

#define  Pi  3.1415 
#define  nss  2 

#define  velocity  300000000 
#define  omega  1885 
#define  dist  500000 
#define  ncc  16 
float  Pm[91]= {23}; 
float  eigenvalues[16]={333}; 
float  eigvects[32][16]={14}; 
float  cova[32][32]={55}; 

_ complex _ float  Part[15]  ={5}; 

_ complex _ float  eigvects2[16][16]  ={5}; 

_ complex _ float  En[  16]  [15]  ={15}; 

_ complex _ float  nos[64][16]  BUFF={5}; 

_ complex _ float  time_delay[16][2]  BUFF  ={2,2, 4, 5, 3, 6, 9, 8, 7, 7, 4, 1,1, 1,1, 5}; 

float  buff21[32][32]  ={45}  ; 
float  buff5[32][32]  ={23}  ; 
float  w[32] ; 
float  rand(); 

//random  number  generator 

int  main() _ attribute _ ((halt, naked)); 

int  randomseed  =0  ; 
int  siggen(); 

//data  generation  function 

_ complex _ float  conj( _ complex _ float  a); 

void  covl(); 

//  Covariance  matrix  computation  parti 
void  cov2(); 

//Covariance  matrix  computation  part2 
void  fourier(); 

void  house(float  x[32],int  k);  //Householder  vector 

function 

void  eigl(int  f); 

void  QR(); 

//QR  decomposition 
void  sym(); 

//Householder  method 
void  QRshift(int  f); 

//  QR  method 
int  count4; 
void  eigfunc(); 
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void  bubbleSort(float  vals[],  int  array_size); 
void  music2(); 
float  Q[32][32]={55}; 
float  R[32][32]  ={55}  ; 
float  H[32][32]={55}; 

_ complex _ float  W[32]  ={ 

#include  "coeff.dat" 


complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 

complex 


float  xntl[64][16]  XM; 
float  xnt2[64][16]  XM; 
float  xnt3[64][16]  XM; 
float  xnt4[64][16]  XM; 
float  xnt5[64][16]  XM; 
float  xnt6[64][16]  XM; 
float  xnt7[64][  16]  XM; 
float  xnt8[64][16]  XM; 
float  xnt9[64][16]  XM; 
float  xntl0[64][16]  XM; 
float  xnt  11  [64]  [16]  XM; 
float  xntl2[64][16]  XM; 
float  xnt  13  [64]  [16]  XM; 
float  xntl4[64][16]  XM; 
float  xntl5[64][16]  XM; 
float  xntl6[64][16]  XM; 
float  xntl7[64][16]  XM; 
float  xntl8[64][16]  XM; 
float  xntl9[64][16]  XM; 
float  xnt20[64][16]  XM; 
float  xnt21[64][16]  XM; 
float  xnt22[64][16]  XM; 
float  xnt23[64][16]  XM; 
float  xnt24[64][16]  XM; 
float  xnt25[64][16]  XM; 
float  xnt26[64][16]  XM; 
float  xnt27[64][16]  XM; 
float  xnt28[64][16]  XM; 
float  xnt29[64][16]  XM; 
float  xnt30[64][16]  XM; 
float  xnt31[64][16]  XM; 
float  xnt32[64][16]  XM; 
float  xnt33[64][16]  XM; 
float  xnt34[64][  16]  XM; 
float  xnt35[64][16]  XM; 
float  xnt36[64][16]  XM; 
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float  xnt37[64][16]  XM; 
float  xnt38[64][16]  XM; 
float  xnt39[64][16]  XM; 
float  xnt40[64][16]  XM; 
float  xnt41[64][16]  XM; 
float  xnt42[64][16]  XM; 
float  xnt43[64][16]  XM; 
float  xnt44[64][16]  XM; 
float  xnt45[64][16]  XM; 
float  xnt46[64][16]  XM; 
float  xnt47[64][  16]  XM; 
float  xnt48[64][16]  XM; 
float  xnt49[64][16]  XM; 
float  xnt50[64][16]  XM; 
float  xnt51[64][16]  XM; 
float  xnt52[64][16]  XM; 
float  xnt53[64][16]  XM; 
float  xnt54[64][16]  XM; 
float  xnt55[64][16]  XM; 
float  xnt56[64][16]  XM; 
float  xnt57[64][  16]  XM; 
float  xnt58[64][16]  XM; 
float  xnt59[64][16]  XM; 
float  xnt60[64][16]  XM; 
float  xnt61[64][16]  XM; 
float  xnt62[64][16]  XM; 
float  xnt63[64][16]  XM; 
float  xnt64[64][16]  XM; 
float  bcov  1  [  1 6]  [  1 6]  XM; 
float  bcov2[16][16]  XM; 
float  bcov3[16][16]  XM; 
float  bcov4[16][16]  XM; 
float  bcov5[16][16]  XM; 
float  bcov6[16][16]  XM; 
float  bcov7[16][16]  XM; 
float  bcov8[16][16]  XM; 
float  bcov9[16][16]  XM; 
float  bcovl0[16][16]  XM; 
float  bcov  1 1  [  1 6]  [  1 6]  XM; 
float  bcovl2[16][16]  XM; 
float  bcov  1 3  [  1 6]  [  1 6]  XM; 
float  bcovl4[16][16]  XM; 
float  bcovl5[16][16]  XM; 
float  bcovl6[16][16]  XM; 
float  bcovl7[16][16]  XM; 
float  bcov  1 8  [  1 6]  [  1 6]  XM; 
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_ complex _ float  bcovl9[16][16]  XM; 

_ complex _ float  bcov20[16][16]  XM; 

_ complex _ float  bcov21[16][16]  XM; 

_ complex _ float  bcov22[16][16]  XM; 

_ complex _ float  bcov23[16][16]  XM; 

_ complex _ float  bcov24[16][16]  XM; 

_ complex _ float  bcov25[16][16]  XM; 

_ complex _ float  bcov26[16][16]  XM; 

_ complex _ float  bcov27[16][16]  XM; 

_ complex _ float  bcov28[16][16]  XM; 

_ complex _ float  bcov29[16][16]  XM; 

_ complex _ float  bcov30[16][16]  XM; 

_ complex _ float  bcov3 1  [  1 6]  [  1 6]  XM; 

_ complex _ float  bcov32[16][16]  XM; 

_ complex _ float  bcov33[16][16]  XM; 

int  main() 

{ 

siggen(); 
fourier(); 
count4  =0; 
covl(); 
cov2(); 
eigl(0); 
sym(); 

QRshift(30); 

eigfunc(); 

bubbleSort(eigenvalues,  32); 

music2(); 

focus(); 

eig2(); 

sym(); 

QRshift(30); 

eigfunc(); 

bubbleSort(eigenvalues,  32); 

music2(); 

return  0; 

} 

_ complex _ float  conj( _ complex _ float  a){ 

_ complex _ float  b; 

_ real _ b= _ real _ a; 

_ imag _ b  =  -1* _ imag _ a; 

return  b;} 
float  rand() 
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{ 

float  test; 

random_seed  =  random_seed  *  25173  +13849; 
//  return  (unsigned  int)(random_seed  ); 
random_seed=(random_seed  %  1024); 
test  =  ((float)(random_seed))/1024; 
return  (test); 

} 

int  siggen()  { 

float  azi[2]={9,50}; 
float  sig_pow[2]={l,l}; 
float  nos_pow=0. 1 ; 
int  count, count2,count3; 
float  temp; 

int  index, index2,index3; 

_ complex float  sig  temp; 

_ complex float  sig[64][2]; 

nos  _pow=sqrt(nos_pow) ; 

for  (count  =0;count<nss;count++){ 

azi[count]=azi[count]  *Pi/ 180; 


} 

nos_pow  =  sqrt(nos_pow); 

for  (count2  =0;count2<ncc;count2++){ 

for  (eount=0;count<nss;count++){ 

temp  =  sin(azi[count])*omega*dist*count2/velocity; 
//temp  =  temp*omega*count2*dist; 


// 


_ real _ time_delay[count2][count]=  cos(temp); 

_ imag _ time_delay  [count2]  [count]=  sin(temp)*- 1 ; 

} 


} 

//  index2=index*64; 

//  index3=  index2+64; 


for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<16;count2++){ 
nos[count]  [count2]=rand(); 
nos[count]  [count2]=nos[count]  [count2]  *nos_pow; 
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} 

} 

for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<2;count2++){ 
sig[count]  [count2]=rand(); 

sig[count]  [count2]=sig  [count]  [count2]  *sig_pow[count2] ; 

} 

} 

for  (count  =  0;count<64;count++)  { 

for  (count2  =  0;count2<16;count2++){ 
sig_temp  =0; 

for  (count3  =  0;count3<2;count3++){ 
sigtemp 

=sig_temp+sig[count][count3]*time_delay[count2][count3]; 

}  ' 

nos[count][count2]  =nos[count][count2]  +sig_temp; 

} 

} 


P32XM_B(nos,xntl ,  1 024); 
return  0; 


} 

void  fourier()  { 

int  countl,count2,count3; 

_ complex _ float  buffi  [64]; 

_ complex _ float  buff2[64]; 

_ complex _ float  tempi [64]; 

XM2P3_B(nos,xntl ,  1 024); 

for  (count  1  =0;countl<16;countl++){ 

for  (count2  =0;count2<64;count2++){ 
buffi  [count2]=nos[count2]  [count  1  ] ; 

//~//  testl[count2]  =buffl[count2]; 

} 

fft64(&W[0],  &buffl [0],  &templ[0],  &buff2[0]); 
for  (count2  =0;count2<64;count2++){ 
nos[count2]  [count  1  ]  =buff2  [count2] ; 

} 

} 


82 


P32XM_B(nos,xntl ,  1 024); 

} 

void  covl(){ 

_ complex _ float  temp[64][16]; 

int  count, count2,count3; 

for  (count=0;count<64;count++){ 

XM2P3_B(nos,xntl ,  1 024); 


for  (count3  =0;count3<16;count3++){ 

temp  [count]  [count3  ]=nos[count4]  [count3  ] ; 

} 


} 

for  (count2  =0;count2<64;count2++){ 

for  (count  =0;count<16;count++){ 

nos[count2]  [count]=temp[count2]  [count] ; 

} 

} 

} 

void  cov2()  { 

int  count, count2,count3; 

_ complex _ float  bcov[16][16]; 

for  (count  =0;count<16;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 
bcov[count]  [count2]=0; 

} 

} 

for  (count=0;count<64;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 

for(count3=0;count3<16;count3++){ 

bcov[count2][count3]=bcov[count2][count3]+(nos[count][count2]*conj(nos[coun 

t][count3])); 
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} 

for  (count  =0;count<16;count++){ 

for(count2=0 ;  count2<  1 6 ;  count2++)  { 

bcov[count][count2]=bcov[count][count2]/64; 
nos[count]  [count2]=bcov[count]  [count2] ; 

} 


} 

P32XM_B(nos,bcovl,64); 

} 

void  eigl(int  f){ 

int  x,y; 

switch(f)  { 

case  0:  XM2P3_B(nos,bcovl,64); 
break; 

case  1 :  XM2P3_B(nos,bcov2,64); 
break; 

case  2:  XM2P3_B(nos,bcov3,64); 
break; 

case  3:  XM2P3_B(nos,bcov4,64); 
break; 

case  4:  XM2P3_B(nos,bcov5,64); 
break; 

case  5:  XM2P3_B(nos,bcov6,64); 
break; 

case  6:  XM2P3_B(nos,bcov7,64); 
break; 

case  7:  XM2P3_B(nos,bcov8,64); 
break; 

case  8:  XM2P3_B(nos,bcov9,64); 
break; 

case  9:  XM2P3_B(nos,bcov  10,64); 


case  10: 

break; 

XM2P3_B(nos,bcovl  1,64); 

case  1 1 : 

break; 

XM2P3_B(nos,bcov  12,64); 

case  12: 

break; 

XM2P3_B(nos,bcov  1 3 ,64); 

case  13: 

break; 

XM2P3_B(nos,bcov  1 4,64); 

case  14: 

break; 

XM2P3_B(nos,bcovl5,64); 

case  15: 

break; 

XM2P3_B(nos,bcov  16,64); 
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case  16: 

break; 

XM2P3_B(nos,bcov  17,64); 

case  17: 

break; 

XM2P3_B(nos,bcovl8,64); 

case  18: 

break; 

XM2P3_B(nos,bcov  19,64); 

case  19: 

break; 

XM2P3_B(nos,bcov20,64); 

case  20: 

break; 

XM2P3_B(nos,bcov2 1 ,64); 

case  21: 

break; 

XM2P3_B(nos,bcov22,64); 

case  22: 

break; 

XM2P3_B(nos,bcov23,64); 

case  23: 

break; 

XM2P3_B(nos,bcov24,64); 

case  24: 

break; 

XM2P3_B(nos,bcov25,64); 

case  25: 

break; 

XM2P3_B(nos,bcov26,64); 

case  26: 

break; 

XM2P3_B(nos,bcov27,64); 

case  27: 

break; 

XM2P3_B(nos,bcov28,64); 

case  28: 

break; 

XM2P3_B(nos,bcov29,64); 

case  29: 

break; 

XM2P3_B(nos,bcov30,64); 

case  30: 

break; 

XM2P3_B(nos,bcov3 1 ,64); 

case  3 1 : 

break; 

XM2P3_B(nos,bcov32,64); 

case  32: 

break; 

XM2P3_B(nos,bcov33,64); 

break; 

} 

//-  //XM2P3_B(m,bcov  1,16); 

for  (x  =0;x<16;x++){ 

for  (y=0;y<16;y++){ 

cova[x][y]  = _ real _ nos[x][y]; 


} 

} 

for  (x  =16;x<32;x++){ 
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for  (y=0;y<16;y++){ 

cova[x][y]  =( imag nos[x-8][y])*(-l); 

} 

} 

for  (x  =0;x<16;x++){ 

for  (y=16;y<32;y++){ 

cova[x][y]  =( imag nos[x][y-8]); 

} 

} 

for  (x  =16;x<32;x++){ 

for  (y=16;y<32;y++){ 

cova[x][y]  = _ real _ nos[x-8][y-8]; 

} 

} 


void  house(float  x[],int  k){ 
float  s, gamma; 
int  i; 

gamma  =  0; 

for  (i  =  0;i<(k-l);i++){ 

w[i]  =  0; 

} 

for  (i  =(k-l);i<32;i++){ 
gamma  +=  x[i]*x[i]; 

} 

gamma  =  sqrt(gamma); 
s  =  x[k-l]*x[k-l]; 
s  =  sqrt(s); 
s  =  gamma+s; 
s  =  2*gamma*s; 
s  =  sqrt(s); 
if  (x[k-l]>=  0){ 

w[k-l]  =  x[k-l]+gamma; 
w[k-l]  =  w[k-l]/s; 

} 
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else  { 

w[k-l]  =  x[k-l] -gamma; 
w[k-l]  =  w[k-l]/s; 

} 

for  (i  =  k;i<32;i++){ 
w[i]  =  x[i]/s; 

} 

} 

void  QR()  { 
float  buff2[32]; 

_ complex _ float  buff[32][32]; 

//float  buff3[32][32]; 
int  i,j,x,y,z; 
for  (i  =  0;i<32;i++){ 

for  (j  =  0;j<32;j++){ 

if  (i  ==j) _ imag _ buflf[i][j]  =1; 

else _ imag _ buff[i][j]  =  0; 

Q[i]D]  = _ imag _ buff[i][j]; 

} 

} 

for  (i  =  0;i<31;i++){ 

for  (j  =0;j<32;j++){ 
buff2[j]  =  cova[j][i]; 

} 

house(buff2,(i+ 1 )) ; 

//gamma  =  0; 
gamma  =  0; 
delta  =0; 
for(x=0;x<n2;x++)  { 
delta+=buff2  [x]  * w[x] ; 

}  for(x=0;x<n2;x++){ 

buff2  [x]  =buff2  [x]  *  delta; 

} 

for  ( x  =0;x<32;x++){ 

for  ( y  =0;y<32;y++){ 


} 

} 

} 

void  sym()  { 


A[x][y]=A[x][y]-buff2[y]; 

Q[x][y]=Q[x][y]-buff2[y]; 


} 
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_ complex _ float  buff[32][32]; 

float  buff2[32]; 
float  buffB  [32]  [32]; 

//~  float  buff4[n2][n2]; 

//~  float  buff6[n2][n2]; 

//~  float  buff7[n2][n2]; 
float  y[32]; 
float  z[32]; 
float  delta; 
int  x,s,i,j,v; 
for  (  x  =0;x<32;x++)  { 

for  (  s  =0;s<32;s++){ 

if  (x  ==s) _ real _ buff[x][s]  =1; 

else _ real _ buff[x][s]  =0; 

buff5[x][s]  = _ real _ buff[x][s]; 

} 

} 

for  ( i  =0;i<30;i++){ 

for  ( j  =0;j<32;j++){ 

buff2[j]  =  cova[j][i]; 
zD']  =  0; 

yD]  =  0; 

} 

delta  =0; 

house(buff2,(i+2)); 

for(x=0;x<32;x++)  { 
delta+=buff2  [x]  * w[x] ; 

} 

for(x=0;x<n2;x++)  { 
buff2  [x]  =buff2  [x]  *  delta; 

} 

for  ( x  =0;x<32;x++){ 

for  ( s  =0;s<32;s++){ 

cova[x]  [s]=cova[x]  [s]  -buff2  [x] ; 
buff5  [x]  [s]=buff5  [x]  [s]  -buff2  [x] ; 

} 

} 

} 

for  (x  =  0;x<32;x++)  { 
for  (  s  =0;s<32;s++){ 
for  (v  =0;  v<32;  v++){ 

buffB  [x]  [s]+=  cova[x]  [v]  *buff5  [v]  [s] ; 

} 

} 

} 

for  (x  =  0;x<n2;x++)  { 
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for  (  s  =0;s<n2;s++){ 

A[x][s]=  cova[x][s]; 

} 

} 

} 

} 

} 

} 

void  QRshift(int  f)  { 

_ complex _ float  buff[32]  [32] ; 

float  buff7[32][32]; 
float  buff2[32]; 
int  x,y,i,z; 

for  (x  =0;x<32;x++){ 

for  (y  =0;y<32;y++){ 

_ imag _ buff[x][y]  =0; 

if  (x  ==y) _ real _ buff[x][y]  =1; 

else _ real _ buff[x][y]  =0; 

} 

} 

for  (i=0;i<f;i++){ 

QR(); 

for  (x  =0;x<32;x++){ 

for  (y  =0;  y<32;  y++)  { 
buff7[x][y]  =0; 

_ imag buff[x][y]  =  0; 

for  (z  =0;  z<32;  z++)  { 

buff7[x][y]  +=  R[x][z]*Q[z][y]; 

_ imag buff[x][y]+= real buff[x]  [z]  *  Q  [z]  [y  ] ; 

} 

//  if  (fabs(buff[x][y])<le-6)  buff[x][y]=0; 

} 

} 

for(  x  =0;x<32;x++){ 

for  (y  =0;  y<32;  y++)  { 

cova[x][y]  =buff7[x][y]; 

_ real buff[x][y]= imag buff[x][y]; 

//  buffo[x][y]= _ imag _ buff[x][y]; 

buff2 1  [x]  [y]= _ imag _ buff[x]  [y] ; 

}  " . " 

} 

} 

int  j; 

for  ( x  =0;x<32;x+=2){ 
j=x/2; 

eigenvalues  [j  ]  =buff7  [x]  [x] ; 
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} 

} 

void  eigfunc()  { 
int  x,s,j,y,z; 
float  test2; 
float  buff23[32][32]; 
for  (  x  =0;x<32;x++)  { 
for  (  s  =0;  s<32;  s++)  { 
buff23[x][s]  =  0; 
for  (  z  =0;  z<32;  z++)  { 

buff23  [x]  [s]+=  buff5  [x]  [z]  *buff2 1  [z]  [s] ; 

} 

//  float  test2  =fabs(buff23[x][s]); 

//  if  (fabs(test2)<le-4)  { 

//  buff23[x][s]  =0;} 

} 

} 

for(  x  =0;x<32;x++){ 

for(  y  =  0;y<32;y+=2){ 
j=y/2; 

eigvects[x]  [j]=buff23  [x]  [y] ; 

} 

} 

} 

void  bubbleSort(float  vals[],  int  array_size) 

{ 

int  i,  j,  k; 
float  tempi; 
float  temp2[32]; 

for  (i  =  (array_size  -  1);  i  >=  0;  i— ) 

{ 

for  (j  =  1 ;  j  <=  i;j++) 

{ 

if  (vals[j-l]  <  vals[j]) 

{ 

tempi  =  vals[j-l]; 
vals[j-l]  =  vals[j]; 
vals[j]  =  tempi; 

for(k  =0;k<32;k++)  { 
temp2  [k] =eigvects  [k]  [j  - 1  ] ; 

} 

for(k  =0;k<32;k++)  { 
eigvects[k]  [j- 1  ]=eigvects[k]  [j] ; 

} 

for(k  =0;k<32;k++)  { 
eigvects  [k]  [j  ]  =temp2  [k] ; 
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} 


} 

} 

} 

for  (i=0;i<16;i++){ 
for  (j=0;j<16;j++){ 

_ real eigvects2  [i]  [j  ]  =  eigvects  [i]  [j  ] ; 

_ imag _ eigvects2[i][j]=  eigvects[i+8][j]; 

} 

} 

} 

void  music2()  { 
int  k,  count; 
int  l,v; 

float  temp,temp2; 
int  kk; 

_ complex _ float  sr; 

_ complex _ float  array[16]={6}; 

_ complex _ float  array2[16]={7}; 

// complex float  Part[15]  ={5}; 

_ complex _ float  Part2[15]  ={5}; 

float  xxx[91]=  {5}; 

_ complex _ float  En2[15][16]  ={5}; 

_ vector int  error; 

for  (1  =  0;l<16;l++){ 

for  (v  =0;v<7;v++){ 

En[l]  [v] =eigvects2  [1]  [v+ 1  ] ; 

En2[v]  [l]=eigvects2[l]  [v+ 1  ] ; 

}  . ' 

} 

for  (k=0;k<=90;k++){ 
xxx[k]=k*Pi; 
xxx[k]  =xxx[k]/180; 
temp  =  sin(xxx[k]); 
temp  =  temp*Pi; 

for  (count  =  0;count<16;count++){ 
temp2  =  temp*count; 

_ real _ array[count]=  cos(temp2); 

_ real _ array2[count]=  cos(temp2); 

_ imag _ array  [count]=  sin(temp2)*-l; 

_ imag _ array2[count]=  sin(temp2)*-l; 

} 

int  s,j; 

for  (s=0;s<7;s++){ 

Part2[s]  =0; 

for  (kk=0;kk<16;kk++){ 
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Part2  [s] +=conj  (En2  [s]  [kk] )  *  array  [kk] ; 

} 

} 

for  <J =0  ;j  <7  ;j ++)  { 

Part[j]  =0; 

for  (kk=0;kk<16;kk++){ 

Part[j]+=conj(array2[kk])*En[kk][j]; 

} 

} 

Pm[k]  =0; 
sr  =0; 

for  fkk=0;kk<7;kk++)  { 

sr+=Part[kk]  *Part2[kk] ; 

} 

Pm[k]  =  ( _ real _ sr  * _ real _ sr)+( _ imag _ sr* _ imag _ sr) 

Pm[k]=sqrt(Pm[k]); 

} 

} 

void  focus()  { 
float  fO; 

float  temp,temp2,temp3; 
inf  t,i,j,i2,j2,k2,12,m2; 

//bl=10.5; 

for  (i  =0;i<nfib;i++)  { 
fbin[i]  =(i+l)*300/16; 

} 

for  (i  =0;i<ncc;i++){ 

for  (j=0;j<ncc;j++){ 

m[i][j]  =0; 

signal[i][j]  =0; 

} 

} 

t  =  (nfib+l)/2; 
ft)  =fbin[t-l]; 
for(i=0;i<nfib;i++)  { 

for(i2=0;i2<nfib;i2++)  { 

for(j  2=0  ;j  2<nfib  ;j  2++)  { 
tl[i2][j2]  =0; 
tpl[i2]Q2]  =0; 
tcl[i2]Q2]  =0; 
tc2[i2][j2]  =0; 
ttl[i2]Q2]  =0; 

} 

} 

omg  =2*Pi*(fO-fbin[i]); 
bl=20; 
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bl  =bl*Pi; 

bl=bl/l  80; 

temp  =sin(bl); 

temp  =sin(bl)/300000000; 

temp  =temp*500000; 

temp  =temp*omg; 

// 

for(k2=0;k2<ncc;k2++)  { 
temp2=  temp*k2; 

_ real _ tl[k2][k2]=  cos(temp2); 

_ imag _ tl  [k2][k2]=  sin(temp2)*- 1 ; 

} 

switch  (i)  { 

case  0: 

XM2P3_B(tpp,bcovl  ,256); 
break; 

case  1: 

XM2P3_B(tpp,bcov2,256); 

break; 

case  2: 

XM2P3_B(tpp,bcov3 ,256); 
break; 

case  3: 

XM2P3_B(tpp,bcov4,256); 

break; 

case  4: 

XM2P3_B(tpp,bcov5 ,256); 
break; 

case  5: 

XM2P3_B(tpp,bcov6,256); 

break; 

case  6: 

XM2P3_B(tpp,bcov7,256); 

break; 

case  7: 

XM2P3_B(tpp,bcov8 ,256); 
break; 

case  8: 

XM2P3_B(tpp,bcov9,256); 

break; 

case  9: 

XM2P3_B(tpp,bcov  10,256); 
break; 
case  10: 

XM2P3_B(tpp,bcovl  1,256); 
break; 
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case  1 1 : 

XM2P3_B(tpp,bcov  12,256); 
break; 
case  12: 

XM2P3_B(tpp,bcov  13,256); 
break; 
case  13: 

XM2P3_B(tpp,bcov  14,256); 
break; 
case  14: 

XM2P3_B(tpp,bcov  15,256); 
break; 
case  15: 

XM2P3_B(tpp,bcov  16,256); 
break; 
case  16: 

XM2P3_B(tpp,bcov  17,256); 
break; 
case  17: 

XM2P3_B(tpp,bcovl  8,256); 
break; 
case  18: 

XM2P3_B(tpp,bcov  19,256); 
break; 
case  19: 

XM2P3_B(tpp,bcov20,256); 
break; 
case  20: 

XM2P3_B(tpp,bcov2 1 ,256); 
break; 
case  21: 

XM2P3_B(tpp,bcov22,256); 
break; 
case  22: 

XM2P3_B(tpp,bcov23 ,256); 

break; 

case  23: 

XM2P3_B(tpp,bcov24,256); 

break; 

case  24: 

XM2P3_B(tpp,bcov25 ,256); 

break; 

case  25: 

XM2P3_B(tpp,bcov26,256); 

break; 

case  26: 
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XM2P3_B(tpp,bcov27,256); 

break; 

case  27: 

XM2P3_B(tpp,bcov28,256); 

break; 

case  28: 

XM2P3_B(tpp,bcov29,256); 

break; 

case  29: 

XM2P3_B(tpp,bcov30,256); 

break; 

case  30: 

XM2P3_B(tpp,bcov3 1 ,256); 

break; 

case  3 1 : 

XM2P3_B(tpp,bcov32,256); 

break; 

case  32: 

XM2P3_B(tpp,bcov3  3,256); 
break; 

} 

for(i2=0;i2<ncc;i2++)  { 

tpl[i2][i2]  =conj(tl[i2][i2]); 
tt  1  [i2]  [i2]=t  1  [i2]  [i2]  *tp  1  [i2]  [i2] ; 

} 

for(  i2  =0;i2<ncc;i2++){ 
for  (  j2  =0;  j2<ncc;  j2++)  { 

tel  [i2][j2]  =  0; 
for  (k2  =0;  k2<ncc;  k2++){ 

tc  1  [i2]  [j2]+=  1 1  [i2]  [k2]  *tpp[k2]  [j2] ; 

} 

} 

} 

for(  i2  =0;i2<ncc;i2++){ 
for  (  j2  =0;  j2<ncc;  j2++)  { 

tc2[i2][j2]  =  0; 
for  (k2  =0;  k2<ncc;  k2++){ 

tc2[i2][j2]+=  tel  [i2][k2]*tpl[k2][j2]; 

} 

signal[i2][j2]  =  signal[i2][j2]+tc2[i2][j2]; 

} 

} 

} 
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